jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1342 from dsbolin/gran_mods 

New generalized granular pair style added
This commit is contained in:
Axel Kohlmeyer 2019-03-29 10:37:07 -04:00 committed by GitHub
parent 93f531441a ab12a7c95b
commit 30929d95e2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
22 changed files with 3665 additions and 698 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_fix.txt
View File

@ -224,7 +224,7 @@ OPT.
"wall/body/polyhedron"_fix_wall_body_polyhedron.html,
"wall/colloid"_fix_wall.html,
"wall/ees"_fix_wall_ees.html,
"wall/gran (o)"_fix_wall_gran.html,
"wall/gran"_fix_wall_gran.html,
"wall/gran/region"_fix_wall_gran_region.html,
"wall/harmonic"_fix_wall.html,
"wall/lj1043"_fix_wall.html,

101
doc/src/fix_wall_gran.txt
View File

@ -7,22 +7,24 @@
:line
fix wall/gran command :h3
fix wall/gran/omp command :h3
[Syntax:]
fix ID group-ID wall/gran fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle args keyword values ... :pre
fix ID group-ID wall/gran fstyle fstyle_params 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
fstyle = style of force interactions between particles and wall :l
possible choices: hooke, hooke/history, hertz/history :pre
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 value between 0.0 and 1.0e4) :l
dampflag = 0 or 1 if tangential damping force is excluded or included :l
possible choices: hooke, hooke/history, hertz/history, granular :pre
fstyle_params = parameters associated with force interaction style :l
For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are:
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below)
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below)
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below)
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)
xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
dampflag = 0 or 1 if tangential damping force is excluded or included :pre
For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre
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
@ -44,7 +46,10 @@ keyword = {wiggle} or {shear} :l
fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0
fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL
fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 :pre
fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0
fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 zplane 0.0 NULL
fix 4 all wall/gran granular jkr 1000.0 50.0 0.3 5.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall zcylinder 15.0 wiggle z 3.0 2.0
fix 5 all wall/gran granular dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall zplane 0.0 NULL :pre
[Description:]
@ -54,31 +59,40 @@ close enough to touch it.
The nature of the wall/particle interactions are determined by the
{fstyle} setting. It can be any of the styles defined by the
"pair_style granular"_pair_gran.html commands. Currently this is
{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). Specifically, delta = radius - r =
overlap of particle with wall, m_eff = mass of particle, and the
effective radius of contact = RiRj/Ri+Rj is just the radius of the
particle.
"pair_style gran/*"_pair_gran.html or the more general "pair_style
granular"_pair_granular.html" commands. Currently the options are
{hooke}, {hooke/history}, or {hertz/history} for the former, and
{granular} with all the possible options of the associated
{pair_coeff} command for the latter. The equation for the force
between the wall and particles touching it is the same as the
corresponding equation on the "pair_style gran/*"_pair_gran.html and
"pair_style_granular"_pair_granular.html doc pages, in the limit of
one of the two particles going to infinite radius and mass (flat
wall). Specifically, delta = radius - r = overlap of particle with
wall, m_eff = mass of particle, and the effective radius of contact =
RiRj/Ri+Rj is set to the radius of the particle.
The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}
have the same meaning and units as those specified with the
"pair_style granular"_pair_gran.html commands. This means a NULL can
be used for either {Kt} or {gamma_t} as described on that page. If a
"pair_style gran/*"_pair_gran.html commands. 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}.
All the model choices for cohesion, tangential friction, rolling
friction and twisting friction supported by the "pair_style
granular"_pair_granular.html through its {pair_coeff} command are also
supported for walls. These are discussed in greater detail on the doc
page for "pair_style granular"_pair_granular.html.
Note that you can choose a different force styles and/or different
values for the 6 wall/particle coefficients than for particle/particle
values for the wall/particle coefficients than for particle/particle
interactions. E.g. if you wish to model the wall as a different
material.
NOTE: As discussed on the doc page for "pair_style
granular"_pair_gran.html, versions of LAMMPS before 9Jan09 used a
gran/*"_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
@ -108,14 +122,14 @@ Optionally, the wall can be moving, if the {wiggle} or {shear}
keywords are appended. Both keywords cannot be used together.
For the {wiggle} keyword, the wall oscillates sinusoidally, similar to
the oscillations of particles which can be specified by the
"fix move"_fix_move.html command. This is useful in packing
simulations of granular particles. The arguments to the {wiggle}
keyword specify a dimension for the motion, as well as it's
{amplitude} and {period}. Note that if the dimension is in the plane
of the wall, this is effectively a shearing motion. If the dimension
is perpendicular to the wall, it is more of a shaking motion. A
{zcylinder} wall can only be wiggled in the z dimension.
the oscillations of particles which can be specified by the "fix
move"_fix_move.html command. This is useful in packing simulations of
granular particles. The arguments to the {wiggle} keyword specify a
dimension for the motion, as well as it's {amplitude} and {period}.
Note that if the dimension is in the plane of the wall, this is
effectively a shearing motion. If the dimension is perpendicular to
the wall, it is more of a shaking motion. A {zcylinder} wall can only
be wiggled in the z dimension.
Each timestep, the position of a wiggled wall in the appropriate {dim}
is set according to this equation:
@ -137,28 +151,6 @@ the clockwise direction for {vshear} > 0 or counter-clockwise for
{vshear} < 0. In this case, {vshear} is the tangential velocity of
the wall at whatever {radius} has been defined.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
[Restart, fix_modify, output, run start/stop, minimize info:]
This fix writes the shear friction state of atoms interacting with the
@ -188,6 +180,7 @@ Any dimension (xyz) that has a granular wall must be non-periodic.
"fix move"_fix_move.html,
"fix wall/gran/region"_fix_wall_gran_region.html,
"pair_style granular"_pair_gran.html
"pair_style gran/*"_pair_gran.html
"pair_style granular"_pair_granular.html
[Default:] none

76
doc/src/fix_wall_gran_region.txt
View File

@ -10,24 +10,30 @@ fix wall/gran/region command :h3
[Syntax:]
fix ID group-ID wall/gran/region fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle regionID :pre
fix ID group-ID wall/gran/region fstyle fstyle_params wallstyle regionID :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
wall/region = style name of this fix command :l
fstyle = style of force interactions between particles and wall :l
possible choices: hooke, hooke/history, hertz/history :pre
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 value between 0.0 and 1.0e4) :l
dampflag = 0 or 1 if tangential damping force is excluded or included :l
possible choices: hooke, hooke/history, hertz/history, granular :pre
fstyle_params = parameters associated with force interaction style :l
For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are:
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below)
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below)
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below)
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)
xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
dampflag = 0 or 1 if tangential damping force is excluded or included :pre
For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre
wallstyle = region (see "fix wall/gran"_fix_wall_gran.html for options for other kinds of walls) :l
region-ID = region whose boundary will act as wall :l,ule
[Examples:]
fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone :pre
fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone
fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 region myBox
fix 4 all wall/gran/region granular jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone
fix 5 all wall/gran/region granular dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall region myCone :pre
[Description:]
@ -42,8 +48,8 @@ Here are snapshots of example models using this command.
Corresponding input scripts can be found in examples/granregion.
Click on the images to see a bigger picture. Movies of these
simulations are "here on the Movies
page"_http://lammps.sandia.gov/movies.html#granregion of the
LAMMPS web site.
page"_http://lammps.sandia.gov/movies.html#granregion of the LAMMPS
web site.
:image(JPG/gran_funnel_small.jpg,JPG/gran_funnel.png)
:image(JPG/gran_mixer_small.jpg,JPG/gran_mixer.png)
@ -123,12 +129,16 @@ to make the two faces differ by epsilon in their position.
The nature of the wall/particle interactions are determined by the
{fstyle} setting. It can be any of the styles defined by the
"pair_style granular"_pair_gran.html commands. Currently this is
{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, but the effective radius is calculated using the radius of the
particle and the radius of curvature of the wall at the contact point.
"pair_style gran/*"_pair_gran.html or the more general "pair_style
granular"_pair_granular.html" commands. Currently the options are
{hooke}, {hooke/history}, or {hertz/history} for the former, and
{granular} with all the possible options of the associated
{pair_coeff} command for the latter. The equation for the force
between the wall and particles touching it is the same as the
corresponding equation on the "pair_style gran/*"_pair_gran.html and
"pair_style_granular"_pair_granular.html doc pages, but the effective
radius is calculated using the radius of the particle and the radius
of curvature of the wall at the contact point.
Specifically, delta = radius - r = overlap of particle with wall,
m_eff = mass of particle, and RiRj/Ri+Rj is the effective radius, with
@ -141,12 +151,18 @@ particle.
The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}
have the same meaning and units as those specified with the
"pair_style granular"_pair_gran.html commands. This means a NULL can
be used for either {Kt} or {gamma_t} as described on that page. If a
"pair_style gran/*"_pair_gran.html commands. 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}.
All the model choices for cohesion, tangential friction, rolling
friction and twisting friction supported by the "pair_style
granular"_pair_granular.html through its {pair_coeff} command are also
supported for walls. These are discussed in greater detail on the doc
page for "pair_style granular"_pair_granular.html.
Note that you can choose a different force styles and/or different
values for the 6 wall/particle coefficients than for particle/particle
interactions. E.g. if you wish to model the wall as a different
@ -154,9 +170,9 @@ material.
[Restart, fix_modify, output, run start/stop, minimize info:]
Similar to "fix wall/gran"_fix_wall_gran.html command, this fix
writes the shear friction state of atoms interacting with the wall to
"binary restart files"_restart.html, so that a simulation can continue
Similar to "fix wall/gran"_fix_wall_gran.html command, this fix writes
the shear friction state of atoms interacting with the wall to "binary
restart files"_restart.html, so that a simulation can continue
correctly if granular potentials with shear "history" effects are
being used. This fix also includes info about a moving region in the
restart file. See the "read_restart"_read_restart.html command for
@ -170,14 +186,14 @@ So you must re-define your region and if it is a moving region, define
its motion attributes in a way that is consistent with the simulation
that wrote the restart file. In particular, if you want to change the
region motion attributes (e.g. its velocity), then you should ensure
the position/orientation of the region at the initial restart
timestep is the same as it was on the timestep the restart file was
written. If this is not possible, you may need to ignore info in the
restart file by defining a new fix wall/gran/region command in your
restart script, e.g. with a different fix ID. Or if you want to keep
the shear history info but discard the region motion information, you
can use the same fix ID for fix wall/gran/region, but assign it a
region with a different region ID.
the position/orientation of the region at the initial restart timestep
is the same as it was on the timestep the restart file was written.
If this is not possible, you may need to ignore info in the restart
file by defining a new fix wall/gran/region command in your restart
script, e.g. with a different fix ID. Or if you want to keep the
shear history info but discard the region motion information, you can
use the same fix ID for fix wall/gran/region, but assign it a region
with a different region ID.
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. No global or per-atom quantities are stored by this fix for

1
doc/src/lammps.book
View File

@ -580,6 +580,7 @@ pair_extep.html
pair_gauss.html
pair_gayberne.html
pair_gran.html
pair_granular.html
pair_gromacs.html
pair_gw.html
pair_ilp_graphene_hbn.html

765
doc/src/pair_granular.txt Normal file
View File

@ -0,0 +1,765 @@
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
pair_style granular command :h3
[Syntax:]
pair_style granular cutoff :pre
cutoff = global cutoff (optional). See discussion below. :l
[Examples:]
pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 :pre
pair_style granular
pair_coeff * * hertz 1000.0 50.0 tangential mindlin NULL 1.0 0.4 :pre
pair_style granular
pair_coeff * * hertz/material 1e8 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji :pre
pair_style granular
pair_coeff 1 1 jkr 1000.0 50.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre
pair_style granular
pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall
pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre
[Description:]
The {granular} styles support a variety of options for the normal,
tangential, rolling and twisting forces resulting from contact between
two granular particles. This expands on the options offered by the
"pair gran/*"_pair_gran.html pair styles. The total computed forces
and torques are the sum of various models selected for the normal,
tangential, rolling and twisting modes of motion.
All model choices and parameters are entered in the
"pair_coeff"_pair_coeff.html command, as described below. Unlike
e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not
global, but can be set to different values for different combinations
of particle types, as determined by the "pair_coeff"_pair_coeff.html
command. If the contact model choice is the same for two particle
types, the mixing for the cross-coefficients can be carried out
automatically. This is shown in the second example, where model
choices are the same for type 1 - type 1 as for type 2 - type2
interactions, but coefficients are different. In this case, the
coefficients for type 2 - type interactions can be determined from
mixing rules discussed below. For additional flexibility,
coefficients as well as model forms can vary between particle types,
as shown in the third example: type 1- type 1 interactions are based
on a Hertzian normal contact model and 2-2 interactions are based on a
DMT cohesive model (see below). In that example, 1-1 and 2-2
interactions have different model forms, in which case mixing of
coefficients cannot be determined, so 1-2 interactions must be
explicitly defined via the {pair_coeff 1 2} command, otherwise an
error would result.
:line
The first required keyword for the {pair_coeff} command is the normal
contact model. Currently supported options for normal contact models
and their required arguments are:
{hooke} : \(k_n\), \(\eta_\{n0\}\) (or \(e\))
{hertz} : \(k_n\), \(\eta_\{n0\}\) (or \(e\))
{hertz/material} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\)
{dmt} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\)
{jkr} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\) :ol
Here, \(k_n\) is spring stiffness (with units that depend on model
choice, see below); \(\eta_\{n0\}\) is a damping prefactor (or, in its
place a coefficient of restitution \(e\), depending on the choice of
damping mode, see below); E is Young's modulus in units of
{force}/{length}^2, i.e. {pressure}; \(\nu\) is Poisson's ratio and
\(\gamma\) is a surface energy density, in units of
{energy}/{length}^2.
For the {hooke} model, the normal, elastic component of force acting
on particle {i} due to contact with particle {j} is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\}
\end\{equation\}
Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle
overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} =
\mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two
particle centers (note the i-j ordering so that \(F_\{ne\}\) is
positive for repulsion), and \(\mathbf\{n\} =
\frac\{\mathbf\{r\}_\{ij\}\}\{\|\mathbf\{r\}_\{ij\}\|\}\). Therefore,
for {hooke}, the units of the spring constant \(k_n\) are
{force}/{distance}, or equivalently {mass}/{time^2}.
For the {hertz} model, the normal component of force is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hertz\} = k_N R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\}
\end\{equation\}
Here, \(R_\{eff\} = \frac\{R_i R_j\}\{R_i + R_j\}\) is the effective
radius, denoted for simplicity as {R} from here on. For {hertz}, the
units of the spring constant \(k_n\) are {force}/{length}^2, or
equivalently {pressure}.
For the {hertz/material} model, the force is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hertz/material\} = \frac\{4\}\{3\} E_\{eff\} R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\}
\end\{equation\}
Here, \(E_\{eff\} = E = \left(\frac\{1-\nu_i^2\}\{E_i\} +
\frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}\) is the effective Young's
modulus, with \(\nu_i, \nu_j \) the Poisson ratios of the particles of
types {i} and {j}. Note that if the elastic modulus and the shear
modulus of the two particles are the same, the {hertz/material} model
is equivalent to the {hertz} model with \(k_N = 4/3 E_\{eff\}\)
The {dmt} model corresponds to the
"(Derjaguin-Muller-Toporov)"_#DMT1975 cohesive model, where the force
is simply Hertz with an additional attractive cohesion term:
\begin\{equation\}
\mathbf\{F\}_\{ne, dmt\} = \left(\frac\{4\}\{3\} E R^\{1/2\}\delta_\{ij\}^\{3/2\} - 4\pi\gamma R\right)\mathbf\{n\}
\end\{equation\}
The {jkr} model is the "(Johnson-Kendall-Roberts)"_#JKR1971 model,
where the force is computed as:
\begin\{equation\}
\label\{eq:force_jkr\}
\mathbf\{F\}_\{ne, jkr\} = \left(\frac\{4Ea^3\}\{3R\} - 2\pi a^2\sqrt\{\frac\{4\gamma E\}\{\pi a\}\}\right)\mathbf\{n\}
\end\{equation\}
Here, {a} is the radius of the contact zone, related to the overlap
\(\delta\) according to:
\begin\{equation\}
\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}
\end\{equation\}
LAMMPS internally inverts the equation above to solve for {a} in terms
of \(\delta\), then solves for the force in the previous
equation. Additionally, note that the JKR model allows for a tensile
force beyond contact (i.e. for \(\delta < 0\)), up to a maximum of
\(3\pi\gamma R\) (also known as the 'pull-off' force). Note that this
is a hysteretic effect, where particles that are not contacting
initially will not experience force until they come into contact
\(\delta \geq 0\); as they move apart and (\(\delta < 0\)), they
experience a tensile force up to \(3\pi\gamma R\), at which point they
lose contact.
:line
In addition, the normal force is augmented by a damping term of the
following general form:
\begin\{equation\}
\mathbf\{F\}_\{n,damp\} = -\eta_n \mathbf\{v\}_\{n,rel\}
\end\{equation\}
Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i)
\cdot \mathbf\{n\}\) is the component of relative velocity along
\(\mathbf\{n\}\).
The optional {damping} keyword to the {pair_coeff} command followed by
a keyword determines the model form of the damping factor \(\eta_n\),
and the interpretation of the \(\eta_\{n0\}\) or \(e\) coefficients
specified as part of the normal contact model settings. The {damping}
keyword and corresponding model form selection may be appended
anywhere in the {pair coeff} command. Note that the choice of damping
model affects both the normal and tangential damping (and depending on
other settings, potentially also the twisting damping). The options
for the damping model currently supported are:
{velocity}
{viscoelastic}
{tsuji} :ol
If the {damping} keyword is not specified, the {viscoelastic} model is
used by default.
For {damping velocity}, the normal damping is simply equal to the
user-specified damping coefficient in the {normal} model:
\begin\{equation\}
\eta_n = \eta_\{n0\}\
\end\{equation\}
Here, \(\gamma_n\) is the damping coefficient specified for the normal
contact model, in units of {mass}/{time},
The {damping viscoelastic} model is based on the viscoelastic
treatment of "(Brilliantov et al)"_#Brill1996, where the normal
damping is given by:
\begin\{equation\}
\eta_n = \eta_\{n0\}\ a m_\{eff\}
\end\{equation\}
Here, \(m_\{eff\} = m_i m_j/(m_i + m_j)\) is the effective mass, {a}
is the contact radius, given by \(a =\sqrt\{R\delta\}\) for all models
except {jkr}, for which it is given implicitly according to \(delta =
a^2/R - 2\sqrt\{\pi \gamma a/E\}\). In this case, \eta_\{n0\}\ is in
units of 1/({time}*{distance}).
The {tsuji} model is based on the work of "(Tsuji et
al)"_#Tsuji1992. Here, the damping coefficient specified as part of
the normal model is interpreted as a restitution coefficient
\(e\). The damping constant \(\eta_n\) is given by:
\begin\{equation\}
\eta_n = \alpha (m_\{eff\}k_n)^\{1/2\}
\end\{equation\}
For normal contact models based on material parameters, \(k_n =
4/3Ea\). The parameter \(\alpha\) is related to the restitution
coefficient {e} according to:
\begin\{equation\}
\alpha = 1.2728-4.2783e+11.087e^2-22.348e^3+27.467e^4-18.022e^5+4.8218e^6
\end\{equation\}
The dimensionless coefficient of restitution \(e\) specified as part
of the normal contact model parameters should be between 0 and 1, but
no error check is performed on this.
The total normal force is computed as the sum of the elastic and
damping components:
\begin\{equation\}
\mathbf\{F\}_n = \mathbf\{F\}_\{ne\} + \mathbf\{F\}_\{n,damp\}
\end\{equation\}
:line
The {pair_coeff} command also requires specification of the tangential
contact model. The required keyword {tangential} is expected, followed
by the model choice and associated parameters. Currently supported
tangential model choices and their expected parameters are as follows:
{linear_nohistory} : \(x_\{\gamma,t\}\), \(\mu_s\)
{linear_history} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\)
{mindlin} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\)
{mindlin_rescale} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) :ol
Here, \(x_\{\gamma,t\}\) is a dimensionless multiplier for the normal
damping \(\eta_n\) that determines the magnitude of the tangential
damping, \(\mu_t\) is the tangential (or sliding) friction
coefficient, and \(k_t\) is the tangential stiffness coefficient.
For {tangential linear_nohistory}, a simple velocity-dependent Coulomb
friction criterion is used, which mimics the behavior of the {pair
gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|\mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
The tangential damping force \(\mathbf\{F\}_\mathrm\{t,damp\}\) is given by:
\begin\{equation\}
\mathbf\{F\}_\mathrm\{t,damp\} = -\eta_t \mathbf\{v\}_\{t,rel\}
\end\{equation\}
The tangential damping prefactor \(\eta_t\) is calculated by scaling
the normal damping \(\eta_n\) (see above):
\begin\{equation\}
\eta_t = -x_\{\gamma,t\} \eta_n
\end\{equation\}
The normal damping prefactor \(\eta_n\) is determined by the choice of
the {damping} keyword, as discussed above. Thus, the {damping}
keyword also affects the tangential damping. The parameter
\(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the
literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
"Tsuji et al"_#Tsuji1992, "Silbert et al"_#Silbert2001). The relative
tangential velocity at the point of contact is given by
\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i +
R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} =
\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r =
\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force
is \(\mathbf\{t\} =
\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\).
The normal force value \(F_\{n0\}\) used to compute the critical force
depends on the form of the contact model. For non-cohesive models
({hertz}, {hertz/material}, {hooke}), it is given by the magnitude of
the normal force:
\begin\{equation\}
F_\{n0\} = \|\mathbf\{F\}_n\|
\end\{equation\}
For cohesive models such as {jkr} and {dmt}, the critical force is
adjusted so that the critical tangential force approaches \(\mu_t
F_\{pulloff\}\), see "Marshall"_#Marshall2009, equation 43, and
"Thornton"_#Thornton1991. For both models, \(F_\{n0\}\) takes the
form:
\begin\{equation\}
F_\{n0\} = \|\mathbf\{F\}_ne + 2 F_\{pulloff\}\|
\end\{equation\}
Where \(F_\{pulloff\} = 3\pi \gamma R \) for {jkr}, and
\(F_\{pulloff\} = 4\pi \gamma R \) for {dmt}.
The remaining tangential options all use accumulated tangential
displacement (i.e. contact history). This is discussed below in the
context of the {linear_history} option, but the same treatment of the
accumulated displacement applies to the other options as well.
For {tangential linear_history}, the tangential force is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t\mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
Here, \(\mathbf\{\xi\}\) is the tangential displacement accumulated
during the entire duration of the contact:
\begin\{equation\}
\mathbf\{\xi\} = \int_\{t0\}^t \mathbf\{v\}_\{t,rel\}(\tau) \mathrm\{d\}\tau
\end\{equation\}
This accumulated tangential displacement must be adjusted to account
for changes in the frame of reference of the contacting pair of
particles during contact. This occurs due to the overall motion of the
contacting particles in a rigid-body-like fashion during the duration
of the contact. There are two modes of motion that are relevant: the
'tumbling' rotation of the contacting pair, which changes the
orientation of the plane in which tangential displacement occurs; and
'spinning' rotation of the contacting pair about the vector connecting
their centers of mass (\(\mathbf\{n\}\)). Corrections due to the
former mode of motion are made by rotating the accumulated
displacement into the plane that is tangential to the contact vector
at each step, or equivalently removing any component of the tangential
displacement that lies along \(\mathbf\{n\}\), and rescaling to
preserve the magnitude. This follows the discussion in
"Luding"_#Luding2008, see equation 17 and relevant discussion in that
work:
\begin\{equation\}
\mathbf\{\xi\} = \left(\mathbf\{\xi'\} - (\mathbf\{n\} \cdot \mathbf\{\xi'\})\mathbf\{n\}\right) \frac\{\|\mathbf\{\xi'\}\|\}\{\|\mathbf\{\xi'\}\| - \mathbf\{n\}\cdot\mathbf\{\xi'\}\}
\label\{eq:rotate_displacements\}
\end\{equation\}
Here, \(\mathbf\{\xi'\}\) is the accumulated displacement prior to the
current time step and \(\mathbf\{\xi\}\) is the corrected
displacement. Corrections to the displacement due to the second mode
of motion described above (rotations about \(\mathbf\{n\}\)) are not
currently implemented, but are expected to be minor for most
simulations.
Furthermore, when the tangential force exceeds the critical force, the
tangential displacement is re-scaled to match the value for the
critical force (see "Luding"_#Luding2008, equation 20 and related
discussion):
\begin\{equation\}
\mathbf\{\xi\} = -\frac\{1\}\{k_t\}\left(\mu_t F_\{n0\}\mathbf\{t\} + \mathbf\{F\}_\{t,damp\}\right)
\end\{equation\}
The tangential force is added to the total normal force (elastic plus
damping) to produce the total force on the particle. The tangential
force also acts at the contact point (defined as the center of the
overlap region) to induce a torque on each particle according to:
\begin\{equation\}
\mathbf\{\tau\}_i = -(R_i - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t
\end\{equation\}
\begin\{equation\}
\mathbf\{\tau\}_j = -(R_j - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t
\end\{equation\}
For {tangential mindlin}, the "Mindlin"_#Mindlin1949 no-slip solution is used, which differs from the {linear_history}
option by an additional factor of {a}, the radius of the contact region. The tangential force is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
Here, {a} is the radius of the contact region, given by \(a = \delta
R\) for all normal contact models, except for {jkr}, where it is given
implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see
discussion above. To match the Mindlin solution, one should set \(k_t
= 8G\), where \(G\) is the shear modulus, related to Young's modulus
\(E\) by \(G = E/(2(1+\nu))\), where \(\nu\) is Poisson's ratio. This
can also be achieved by specifying {NULL} for \(k_t\), in which case a
normal contact model that specifies material parameters \(E\) and
\(\nu\) is required (e.g. {hertz/material}, {dmt} or {jkr}). In this
case, mixing of the shear modulus for different particle types {i} and
{j} is done according to:
\begin\{equation\}
1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j
\end\{equation\}
The {mindlin_rescale} option uses the same form as {mindlin}, but the
magnitude of the tangential displacement is re-scaled as the contact
unloads, i.e. if \(a < a_\{t_\{n-1\}\}\):
\begin\{equation\}
\mathbf\{\xi\} = \mathbf\{\xi_\{t_\{n-1\}\}\} \frac\{a\}\{a_\{t_\{n-1\}\}\}
\end\{equation\}
Here, \(t_\{n-1\}\) indicates the value at the previous time
step. This rescaling accounts for the fact that a decrease in the
contact area upon unloading leads to the contact being unable to
support the previous tangential loading, and spurious energy is
created without the rescaling above ("Walton"_#WaltonPC ). See also
discussion in "Thornton et al, 2013"_#Thornton2013 , particularly
equation 18(b) of that work and associated discussion.
:line
The optional {rolling} keyword enables rolling friction, which resists
pure rolling motion of particles. The options currently supported are:
{none}
{sds} : \(k_\{roll\}\), \(\gamma_\{roll\}\), \(\mu_\{roll\}\) :ol
If the {rolling} keyword is not specified, the model defaults to {none}.
For {rolling sds}, rolling friction is computed via a
spring-dashpot-slider, using a 'pseudo-force' formulation, as detailed
by "Luding"_#Luding2008. Unlike the formulation in
"Marshall"_#Marshall2009, this allows for the required adjustment of
rolling displacement due to changes in the frame of reference of the
contacting pair. The rolling pseudo-force is computed analogously to
the tangential force:
\begin\{equation\}
\mathbf\{F\}_\{roll,0\} = k_\{roll\} \mathbf\{\xi\}_\{roll\} - \gamma_\{roll\} \mathbf\{v\}_\{roll\}
\end\{equation\}
Here, \(\mathbf\{v\}_\{roll\} = -R(\mathbf\{\Omega\}_i -
\mathbf\{\Omega\}_j) \times \mathbf\{n\}\) is the relative rolling
velocity, as given in "Wang et al"_#Wang2015 and
"Luding"_#Luding2008. This differs from the expressions given by "Kuhn
and Bagi"_#Kuhn2004 and used in "Marshall"_#Marshall2009; see "Wang et
al"_#Wang2015 for details. The rolling displacement is given by:
\begin\{equation\}
\mathbf\{\xi\}_\{roll\} = \int_\{t_0\}^t \mathbf\{v\}_\{roll\} (\tau) \mathrm\{d\} \tau
\end\{equation\}
A Coulomb friction criterion truncates the rolling pseudo-force if it
exceeds a critical value:
\begin\{equation\}
\mathbf\{F\}_\{roll\} = min(\mu_\{roll\} F_\{n,0\}, \|\mathbf\{F\}_\{roll,0\}\|)\mathbf\{k\}
\end\{equation\}
Here, \(\mathbf\{k\} =
\mathbf\{v\}_\{roll\}/\|\mathbf\{v\}_\{roll\}\|\) is the direction of
the pseudo-force. As with tangential displacement, the rolling
displacement is rescaled when the critical force is exceeded, so that
the spring length corresponds the critical force. Additionally, the
displacement is adjusted to account for rotations of the frame of
reference of the two contacting particles in a manner analogous to the
tangential displacement.
The rolling pseudo-force does not contribute to the total force on
either particle (hence 'pseudo'), but acts only to induce an equal and
opposite torque on each particle, according to:
\begin\{equation\}
\tau_\{roll,i\} = R_\{eff\} \mathbf\{n\} \times \mathbf\{F\}_\{roll\}
\end\{equation\}
\begin\{equation\}
\tau_\{roll,j\} = -\tau_\{roll,i\}
\end\{equation\}
:line
The optional {twisting} keyword enables twisting friction, which
resists rotation of two contacting particles about the vector
\(\mathbf\{n\}\) that connects their centers. The options currently
supported are:
{none}
{sds} : \(k_\{twist\}\), \(\gamma_\{twist\}\), \(\mu_\{twist\}\)
{marshall} :ol
If the {twisting} keyword is not specified, the model defaults to {none}.
For both {twisting sds} and {twisting marshall}, a history-dependent
spring-dashpot-slider is used to compute the twisting torque. Because
twisting displacement is a scalar, there is no need to adjust for
changes in the frame of reference due to rotations of the particle
pair. The formulation in "Marshall"_#Marshall2009 therefore provides
the most straightforward treatment:
\begin\{equation\}
\tau_\{twist,0\} = -k_\{twist\}\xi_\{twist\} - \gamma_\{twist\}\Omega_\{twist\}
\end\{equation\}
Here \(\xi_\{twist\} = \int_\{t_0\}^t \Omega_\{twist\} (\tau)
\mathrm\{d\}\tau\) is the twisting angular displacement, and
\(\Omega_\{twist\} = (\mathbf\{\Omega\}_i - \mathbf\{\Omega\}_j) \cdot
\mathbf\{n\}\) is the relative twisting angular velocity. The torque
is then truncated according to:
\begin\{equation\}
\tau_\{twist\} = min(\mu_\{twist\} F_\{n,0\}, \tau_\{twist,0\})
\end\{equation\}
Similar to the sliding and rolling displacement, the angular
displacement is rescaled so that it corresponds to the critical value
if the twisting torque exceeds this critical value:
\begin\{equation\}
\xi_\{twist\} = \frac\{1\}\{k_\{twist\}\} (\mu_\{twist\} F_\{n,0\}sgn(\Omega_\{twist\}) - \gamma_\{twist\}\Omega_\{twist\})
\end\{equation\}
For {twisting sds}, the coefficients \(k_\{twist\}, \gamma_\{twist\}\)
and \(\mu_\{twist\}\) are simply the user input parameters that follow
the {twisting sds} keywords in the {pair_coeff} command.
For {twisting_marshall}, the coefficients are expressed in terms of
sliding friction coefficients, as discussed in
"Marshall"_#Marshall2009 (see equations 32 and 33 of that work):
\begin\{equation\}
k_\{twist\} = 0.5k_ta^2
\end\{equation\}
\begin\{equation\}
\eta_\{twist\} = 0.5\eta_ta^2
\end\{equation\}
\begin\{equation\}
\mu_\{twist\} = \frac\{2\}\{3\}a\mu_t
\end\{equation\}
Finally, the twisting torque on each particle is given by:
\begin\{equation\}
\mathbf\{\tau\}_\{twist,i\} = \tau_\{twist\}\mathbf\{n\}
\end\{equation\}
\begin\{equation\}
\mathbf\{\tau\}_\{twist,j\} = -\mathbf\{\tau\}_\{twist,i\}
\end\{equation\}
:line
LAMMPS automatically sets pairwise cutoff values for {pair_style
granular} based on particle radii (and in the case of {jkr} pull-off
distances). In the vast majority of situations, this is adequate.
However, a cutoff value can optionally be appended to the {pair_style
granular} command to specify a global cutoff (i.e. a cutoff for all
atom types). Additionally, the optional {cutoff} keyword can be passed
to the {pair_coeff} command, followed by a cutoff value. This will
set a pairwise cutoff for the atom types in the {pair_coeff} command.
These options may be useful in some rare cases where the automatic
cutoff determination is not sufficient, e.g. if particle diameters
are being modified via the {fix adapt} command. In that case, the
global cutoff specified as part of the {pair_style granular} command
is applied to all atom types, unless it is overridden for a given atom
type combination by the {cutoff} value specified in the {pair coeff}
command. If {cutoff} is only specified in the {pair coeff} command
and no global cutoff is appended to the {pair_style granular} command,
then LAMMPS will use that cutoff for the specified atom type
combination, and automatically set pairwise cutoffs for the remaining
atom types.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
The "pair_modify"_pair_modify.html mix, shift, table, and tail options
are not relevant for granular pair styles.
Mixing of coefficients is carried out using geometric averaging for
most quantities, e.g. if friction coefficient for type 1-type 1
interactions is set to \(\mu_1\), and friction coefficient for type
2-type 2 interactions is set to \(\mu_2\), the friction coefficient
for type1-type2 interactions is computed as \(\sqrt\{\mu_1\mu_2\}\)
(unless explicitly specified to a different value by a {pair_coeff 1 2
...} command. The exception to this is elastic modulus, only
applicable to {hertz/material}, {dmt} and {jkr} normal contact
models. In that case, the effective elastic modulus is computed as:
\begin\{equation\}
E_\{eff,ij\} = \left(\frac\{1-\nu_i^2\}\{E_i\} + \frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}
\end\{equation\}
If the {i-j} coefficients \(E_\{ij\}\) and \(\nu_\{ij\}\) are
explicitly specified, the effective modulus is computed as:
\begin\{equation\}
E_\{eff,ij\} = \left(\frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\} + \frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\}\right)^\{-1\}
\end\{equation\}
or
\begin\{equation\}
E_\{eff,ij\} = \frac\{E_\{ij\}\}\{2(1-\nu_\{ij\})\}
\end\{equation\}
These pair styles write their information to "binary restart
files"_restart.html, so a pair_style command does not need to be
specified in an input script that reads a restart file.
These pair styles can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. They do not support the
{inner}, {middle}, {outer} keywords.
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 10 extra pairwise quantities. The first 3 are the
components of the tangential force between particles I and J, acting
on particle I. The 4th 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 last 3 (10-12) are the components of the vector connecting
the centers of the two particles (x_I - x_J).
These extra quantities can be accessed by the "compute
pair/local"_compute_pair_local.html command, as {p1}, {p2}, ...,
{p12}.
:line
[Restrictions:]
All the granular pair styles are part of the GRANULAR package. It is
only enabled if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
These pair styles require that atoms store torque and angular velocity
(omega) as defined by the "atom_style"_atom_style.html. They also
require a per-particle radius is stored. The {sphere} atom style does
all of this.
This pair style requires you to use the "comm_modify vel
yes"_comm_modify.html command so that velocities are stored by ghost
atoms.
These pair styles will not restart exactly when using the
"read_restart"_read_restart.html command, though they should provide
statistically similar results. This is because the forces they
compute depend on atom velocities. See the
"read_restart"_read_restart.html command for more details.
[Related commands:]
"pair_coeff"_pair_coeff.html
"pair gran/*"_pair_gran.html
[Default:]
For the {pair_coeff} settings: {damping viscoelastic}, {rolling none},
{twisting none}.
[References:]
:link(Brill1996)
[(Brilliantov et al, 1996)] Brilliantov, N. V., Spahn, F., Hertzsch,
J. M., & Poschel, T. (1996). Model for collisions in granular
gases. Physical review E, 53(5), 5382.
:link(Tsuji1992)
[(Tsuji et al, 1992)] Tsuji, Y., Tanaka, T., & Ishida,
T. (1992). Lagrangian numerical simulation of plug flow of
cohesionless particles in a horizontal pipe. Powder technology, 71(3),
239-250.
:link(JKR1971)
[(Johnson et al, 1971)] Johnson, K. L., Kendall, K., & Roberts,
A. D. (1971). Surface energy and the contact of elastic
solids. Proc. R. Soc. Lond. A, 324(1558), 301-313.
:link(DMT1975)
[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.
:link(Luding2008)
[(Luding, 2008)] Luding, S. (2008). Cohesive, frictional powders:
contact models for tension. Granular matter, 10(4), 235.
:link(Marshall2009)
[(Marshall, 2009)] Marshall, J. S. (2009). Discrete-element modeling
of particulate aerosol flows. Journal of Computational Physics,
228(5), 1541-1561.
:link(Silbert2001)
[(Silbert, 2001)] Silbert, L. E., Ertas, D., Grest, G. S., Halsey,
T. C., Levine, D., & Plimpton, S. J. (2001). Granular flow down an
inclined plane: Bagnold scaling and rheology. Physical Review E,
64(5), 051302.
:link(Kuhn2004)
[(Kuhn and Bagi, 2005)] Kuhn, M. R., & Bagi, K. (2004). Contact
rolling and deformation in granular media. International journal of
solids and structures, 41(21), 5793-5820.
:link(Wang2015)
[(Wang et al, 2015)] Wang, Y., Alonso-Marroquin, F., & Guo,
W. W. (2015). Rolling and sliding in 3-D discrete element
models. Particuology, 23, 49-55.
:link(Thornton1991)
[(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the
presence of adhesion. J. Phys. D: Appl. Phys. 24 1942
:link(Mindlin1949)
[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies
in contact. J. Appl. Mech., ASME 16, 259-268.
:link(Thornton2013)
[(Thornton et al, 2013)] Thornton, C., Cummins, S. J., & Cleary,
P. W. (2013). An investigation of the comparative behaviour of
alternative contact force models during inelastic collisions. Powder
Technology, 233, 30-46.
:link(WaltonPC)
[(Otis R. Walton)] Walton, O.R., Personal Communication

1
doc/src/pairs.txt
View File

@ -42,6 +42,7 @@ Pair Styles :h1
pair_gauss
pair_gayberne
pair_gran
pair_granular
pair_gromacs
pair_gw
pair_hbond_dreiding

19
doc/utils/sphinx-config/false_positives.txt
View File

@ -156,6 +156,8 @@ ba
Babadi
backcolor
Baczewski
Bagi
Bagnold
Bal
balancer
Balankura
@ -347,6 +349,7 @@ Cij
cis
civ
clearstore
Cleary
Clebsch
clemson
Clermont
@ -373,6 +376,7 @@ Coeff
CoefficientN
coeffs
Coeffs
cohesionless
Coker
Colberg
coleman
@ -450,6 +454,7 @@ cuda
Cuda
CUDA
CuH
Cummins
Curk
customIDs
cutbond
@ -493,6 +498,7 @@ darkturquoise
darkviolet
Das
Dasgupta
dashpot
dat
datafile
datums
@ -530,6 +536,7 @@ Dequidt
der
derekt
Derjagin
Derjaguin
Derlet
Deserno
Destree
@ -1081,6 +1088,7 @@ Hyoungki
hyperdynamics
hyperradius
hyperspherical
hysteretic
Ibanez
ibar
ibm
@ -1140,6 +1148,7 @@ interconvert
interial
interlayer
intermolecular
Interparticle
interstitials
Intr
intra
@ -1158,6 +1167,7 @@ IPython
Isele
isenthalpic
ish
Ishida
iso
isodemic
isoenergetic
@ -1453,6 +1463,7 @@ logfile
logfreq
logicals
Lomdahl
Lond
lookups
Lookups
LoopVar
@ -1468,6 +1479,7 @@ lsfftw
ltbbmalloc
lubricateU
lucy
Luding
Lussetti
Lustig
lwsock
@ -1506,6 +1518,7 @@ manybody
MANYBODY
Maras
Marrink
Marroquin
Marsaglia
Marseille
Martyna
@ -1517,6 +1530,7 @@ masstotal
Masuhiro
Matchett
Materias
mathbf
matlab
matplotlib
Mattox
@ -1605,6 +1619,7 @@ Mie
Mikami
Militzer
Minary
Mindlin
mincap
mingw
minima
@ -2290,6 +2305,7 @@ rg
Rg
Rhaphson
rheological
rheology
rhodo
Rhodo
rhodopsin
@ -2606,6 +2622,7 @@ Tait
taitwater
Tajkhorshid
Tamaskovics
Tanaka
tanh
Tartakovsky
taskset
@ -2693,6 +2710,7 @@ tokyo
tol
toolchain
topologies
Toporov
Torder
torsions
Tosi
@ -2737,6 +2755,7 @@ Tsrd
Tstart
tstat
Tstop
Tsuji
Tsuzuki
tt
Tt

1042
src/GRANULAR/fix_wall_gran.cpp
View File

File diff suppressed because it is too large Load Diff

58
src/GRANULAR/fix_wall_gran.h
View File

@ -46,42 +46,70 @@ class FixWallGran : public Fix {
virtual int maxsize_restart();
void reset_dt();
void hooke(double, double, double, double, double *,
double *, double *, double *, double *, double, double);
void hooke(double, double, double, double, double *, double *,
double *, double *, double *, double, double, double*);
void hooke_history(double, double, double, double, double *,
double *, double *, double *, double *, double, double,
double *);
void hertz_history(double, double, double, double, double *, double,
double *, double *, double *, double *, double, double,
double *);
void bonded_history(double, double, double, double, double *, double,
double *, double *, double *, double *, double, double,
double *);
double *, double *, double *, double *, double,
double, double *, double *);
void hertz_history(double, double, double, double, double *,
double, double *, double *, double *, double *,
double, double, double *, double *);
void granular(double, double, double, double, double *, double,
double *, double *, double *, double *, double,
double, double *, double *);
double pulloff_distance(double);
protected:
int wallstyle,wiggle,wshear,axis;
int pairstyle,nlevels_respa;
bigint time_origin;
double kn,kt,gamman,gammat,xmu;
double E,G,SurfEnergy;
// for granular model choices
int normal_model, damping_model;
int tangential_model, roll_model, twist_model;
// history flags
int normal_history, tangential_history, roll_history, twist_history;
// indices of history entries
int normal_history_index;
int tangential_history_index;
int roll_history_index;
int twist_history_index;
// material coefficients
double Emod, poiss, Gmod;
// contact model coefficients
double normal_coeffs[4];
double tangential_coeffs[3];
double roll_coeffs[3];
double twist_coeffs[3];
double lo,hi,cylradius;
double amplitude,period,omega,vshear;
double dt;
char *idregion;
int history; // if particle/wall interaction stores history
int shearupdate; // flag for whether shear history is updated
int sheardim; // # of shear history values per contact
int use_history; // if particle/wall interaction stores history
int history_update; // flag for whether shear history is updated
int size_history; // # of shear history values per contact
// shear history for single contact per particle
double **shearone;
double **history_one;
// rigid body masses for use in granular interactions
class Fix *fix_rigid; // ptr to rigid body fix, NULL if none
double *mass_rigid; // rigid mass for owned+ghost atoms
int nmax; // allocated size of mass_rigid
// store particle interactions
int store;
};
}

186
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -39,15 +39,17 @@ using namespace MathConst;
// same as FixWallGran
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY};
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,GRANULAR};
enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR};
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), ncontact(NULL),
walls(NULL), shearmany(NULL), c2r(NULL)
FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL),
ncontact(NULL),
walls(NULL), history_many(NULL), c2r(NULL)
{
restart_global = 1;
motion_resetflag = 0;
@ -66,17 +68,17 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
// re-allocate atom-based arrays with nshear
// do not register with Atom class, since parent class did that
memory->destroy(shearone);
shearone = NULL;
memory->destroy(history_one);
history_one = NULL;
ncontact = NULL;
walls = NULL;
shearmany = NULL;
history_many = NULL;
grow_arrays(atom->nmax);
// initialize shear history as if particle is not touching region
if (history) {
if (use_history) {
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
ncontact[i] = 0;
@ -92,7 +94,7 @@ FixWallGranRegion::~FixWallGranRegion()
memory->destroy(ncontact);
memory->destroy(walls);
memory->destroy(shearmany);
memory->destroy(history_many);
}
/* ---------------------------------------------------------------------- */
@ -138,8 +140,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
// do not update shear history during setup
shearupdate = 1;
if (update->setupflag) shearupdate = 0;
history_update = 1;
if (update->setupflag) history_update = 0;
// if just reneighbored:
// update rigid body masses for owned atoms if using FixRigid
@ -188,7 +190,13 @@ void FixWallGranRegion::post_force(int /*vflag*/)
if (mask[i] & groupbit) {
if (!region->match(x[i][0],x[i][1],x[i][2])) continue;
nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]);
if (pairstyle == GRANULAR && normal_model == JKR){
nc = region->surface(x[i][0],x[i][1],x[i][2],
radius[i]+pulloff_distance(radius[i]));
}
else{
nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]);
}
if (nc > tmax)
error->one(FLERR,"Too many wall/gran/region contacts for one particle");
@ -198,7 +206,7 @@ void FixWallGranRegion::post_force(int /*vflag*/)
// also set c2r[] = indices into region->contact[] for each of N contacts
// process zero or one contact here, otherwise invoke update_contacts()
if (history) {
if (use_history) {
if (nc == 0) {
ncontact[i] = 0;
continue;
@ -209,15 +217,14 @@ void FixWallGranRegion::post_force(int /*vflag*/)
if (ncontact[i] == 0) {
ncontact[i] = 1;
walls[i][0] = iwall;
for (m = 0; m < sheardim; m++)
shearmany[i][0][m] = 0.0;
for (m = 0; m < size_history; m++)
history_many[i][0][m] = 0.0;
} else if (ncontact[i] > 1 || iwall != walls[i][0])
update_contacts(i,nc);
} else update_contacts(i,nc);
}
// process current contacts
for (int ic = 0; ic < nc; ic++) {
// rsq = squared contact distance
@ -225,36 +232,57 @@ void FixWallGranRegion::post_force(int /*vflag*/)
rsq = region->contact[ic].r*region->contact[ic].r;
if (pairstyle == GRANULAR && normal_model == JKR){
if (history_many[i][c2r[ic]][0] == 0.0 && rsq > radius[i]*radius[i]){
for (m = 0; m < size_history; m++)
history_many[i][0][m] = 0.0;
continue;
}
}
dx = region->contact[ic].delx;
dy = region->contact[ic].dely;
dz = region->contact[ic].delz;
if (regiondynamic) region->velocity_contact(vwall, x[i], ic);
// meff = effective mass of sphere
// if I is part of rigid body, use body mass
meff = rmass[i];
if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i];
// store contact info
if (peratom_flag){
array_atom[i][0] = (double)atom->tag[i];
array_atom[i][4] = x[i][0] - dx;
array_atom[i][5] = x[i][1] - dy;
array_atom[i][6] = x[i][2] - dz;
array_atom[i][7] = radius[i];
}
// invoke sphere/wall interaction
double *contact;
if (peratom_flag)
contact = array_atom[i];
else
contact = NULL;
if (pairstyle == HOOKE)
hooke(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff);
omega[i],torque[i],radius[i],meff, contact);
else if (pairstyle == HOOKE_HISTORY)
hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,
shearmany[i][c2r[ic]]);
omega[i],torque[i],radius[i],meff,
history_many[i][c2r[ic]], contact);
else if (pairstyle == HERTZ_HISTORY)
hertz_history(rsq,dx,dy,dz,vwall,region->contact[ic].radius,
v[i],f[i],omega[i],torque[i],
radius[i],meff,shearmany[i][c2r[ic]]);
else if (pairstyle == BONDED_HISTORY)
bonded_history(rsq,dx,dy,dz,vwall,region->contact[ic].radius,
v[i],f[i],omega[i],torque[i],
radius[i],meff,shearmany[i][c2r[ic]]);
v[i],f[i],omega[i],torque[i],
radius[i],meff,history_many[i][c2r[ic]], contact);
else if (pairstyle == GRANULAR)
granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius,
v[i],f[i],omega[i],torque[i],
radius[i],meff,history_many[i][c2r[ic]],contact);
}
}
}
@ -282,8 +310,8 @@ void FixWallGranRegion::update_contacts(int i, int nc)
if (region->contact[m].iwall == walls[i][iold]) break;
if (m >= nc) {
ilast = ncontact[i]-1;
for (j = 0; j < sheardim; j++)
shearmany[i][iold][j] = shearmany[i][ilast][j];
for (j = 0; j < size_history; j++)
history_many[i][iold][j] = history_many[i][ilast][j];
walls[i][iold] = walls[i][ilast];
ncontact[i]--;
} else iold++;
@ -305,8 +333,8 @@ void FixWallGranRegion::update_contacts(int i, int nc)
iadd = ncontact[i];
c2r[iadd] = inew;
for (j = 0; j < sheardim; j++)
shearmany[i][iadd][j] = 0.0;
for (j = 0; j < size_history; j++)
history_many[i][iadd][j] = 0.0;
walls[i][iadd] = iwall;
ncontact[i]++;
}
@ -321,10 +349,10 @@ double FixWallGranRegion::memory_usage()
{
int nmax = atom->nmax;
double bytes = 0.0;
if (history) { // shear history
if (use_history) { // shear history
bytes += nmax * sizeof(int); // ncontact
bytes += nmax*tmax * sizeof(int); // walls
bytes += nmax*tmax*sheardim * sizeof(double); // shearmany
bytes += nmax*tmax*size_history * sizeof(double); // history_many
}
if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid
return bytes;
@ -336,11 +364,14 @@ double FixWallGranRegion::memory_usage()
void FixWallGranRegion::grow_arrays(int nmax)
{
if (history) {
if (use_history) {
memory->grow(ncontact,nmax,"fix_wall_gran:ncontact");
memory->grow(walls,nmax,tmax,"fix_wall_gran:walls");
memory->grow(shearmany,nmax,tmax,sheardim,"fix_wall_gran:shearmany");
memory->grow(history_many,nmax,tmax,size_history,
"fix_wall_gran:history_many");
}
if (peratom_flag)
memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");
}
/* ----------------------------------------------------------------------
@ -351,16 +382,20 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/)
{
int m,n,iwall;
if (!history) return;
n = ncontact[i];
for (iwall = 0; iwall < n; iwall++) {
walls[j][iwall] = walls[i][iwall];
for (m = 0; m < sheardim; m++)
shearmany[j][iwall][m] = shearmany[i][iwall][m];
if (use_history){
n = ncontact[i];
for (iwall = 0; iwall < n; iwall++) {
walls[j][iwall] = walls[i][iwall];
for (m = 0; m < size_history; m++)
history_many[j][iwall][m] = history_many[i][iwall][m];
}
ncontact[j] = ncontact[i];
}
if (peratom_flag){
for (int m = 0; m < size_peratom_cols; m++)
array_atom[j][m] = array_atom[i][m];
}
ncontact[j] = ncontact[i];
}
/* ----------------------------------------------------------------------
@ -369,8 +404,12 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/)
void FixWallGranRegion::set_arrays(int i)
{
if (!history) return;
ncontact[i] = 0;
if (use_history)
ncontact[i] = 0;
if (peratom_flag){
for (int m = 0; m < size_peratom_cols; m++)
array_atom[i][m] = 0;
}
}
/* ----------------------------------------------------------------------
@ -381,16 +420,19 @@ int FixWallGranRegion::pack_exchange(int i, double *buf)
{
int m;
if (!history) return 0;
int n = 0;
int count = ncontact[i];
buf[n++] = ubuf(count).d;
for (int iwall = 0; iwall < count; iwall++) {
buf[n++] = ubuf(walls[i][iwall]).d;
for (m = 0; m < sheardim; m++)
buf[n++] = shearmany[i][iwall][m];
if (use_history){
int count = ncontact[i];
buf[n++] = ubuf(count).d;
for (int iwall = 0; iwall < count; iwall++) {
buf[n++] = ubuf(walls[i][iwall]).d;
for (m = 0; m < size_history; m++)
buf[n++] = history_many[i][iwall][m];
}
}
if (peratom_flag){
for (int m = 0; m < size_peratom_cols; m++)
buf[n++] = array_atom[i][m];
}
return n;
@ -404,15 +446,19 @@ int FixWallGranRegion::unpack_exchange(int nlocal, double *buf)
{
int m;
if (!history) return 0;
int n = 0;
int count = ncontact[nlocal] = (int) ubuf(buf[n++]).i;
for (int iwall = 0; iwall < count; iwall++) {
walls[nlocal][iwall] = (int) ubuf(buf[n++]).i;
for (m = 0; m < sheardim; m++)
shearmany[nlocal][iwall][m] = buf[n++];
if (use_history){
int count = ncontact[nlocal] = (int) ubuf(buf[n++]).i;
for (int iwall = 0; iwall < count; iwall++) {
walls[nlocal][iwall] = (int) ubuf(buf[n++]).i;
for (m = 0; m < size_history; m++)
history_many[nlocal][iwall][m] = buf[n++];
}
}
if (peratom_flag){
for (int m = 0; m < size_peratom_cols; m++)
array_atom[nlocal][m] = buf[n++];
}
return n;
@ -426,7 +472,7 @@ int FixWallGranRegion::pack_restart(int i, double *buf)
{
int m;
if (!history) return 0;
if (!use_history) return 0;
int n = 1;
int count = ncontact[i];
@ -434,8 +480,8 @@ int FixWallGranRegion::pack_restart(int i, double *buf)
buf[n++] = ubuf(count).d;
for (int iwall = 0; iwall < count; iwall++) {
buf[n++] = ubuf(walls[i][iwall]).d;
for (m = 0; m < sheardim; m++)
buf[n++] = shearmany[i][iwall][m];
for (m = 0; m < size_history; m++)
buf[n++] = history_many[i][iwall][m];
}
buf[0] = n;
return n;
@ -449,7 +495,7 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
{
int k;
if (!history) return;
if (!use_history) return;
double **extra = atom->extra;
@ -462,8 +508,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
int count = ncontact[nlocal] = (int) ubuf(extra[nlocal][m++]).i;
for (int iwall = 0; iwall < count; iwall++) {
walls[nlocal][iwall] = (int) ubuf(extra[nlocal][m++]).i;
for (k = 0; k < sheardim; k++)
shearmany[nlocal][iwall][k] = extra[nlocal][m++];
for (k = 0; k < size_history; k++)
history_many[nlocal][iwall][k] = extra[nlocal][m++];
}
}
@ -473,8 +519,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
int FixWallGranRegion::maxsize_restart()
{
if (!history) return 0;
return 2 + tmax*(sheardim+1);
if (!use_history) return 0;
return 2 + tmax*(size_history+1);
}
/* ----------------------------------------------------------------------
@ -483,8 +529,8 @@ int FixWallGranRegion::maxsize_restart()
int FixWallGranRegion::size_restart(int nlocal)
{
if (!history) return 0;
return 2 + ncontact[nlocal]*(sheardim+1);
if (!use_history) return 0;
return 2 + ncontact[nlocal]*(size_history+1);
}
/* ----------------------------------------------------------------------

2
src/GRANULAR/fix_wall_gran_region.h
View File

@ -54,7 +54,7 @@ class FixWallGranRegion : public FixWallGran {
int tmax; // max # of region walls one particle can touch
int *ncontact; // # of shear contacts per particle
int **walls; // which wall each contact is with
double ***shearmany; // shear history per particle per contact
double ***history_many; // history per particle per contact
int *c2r; // contact to region mapping
// c2r[i] = index of Ith contact in
// region-contact[] list of contacts

7
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -44,6 +44,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
single_enable = 1;
no_virial_fdotr_compute = 1;
history = 1;
size_history = 3;
fix_history = NULL;
single_extra = 10;
@ -57,6 +58,10 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
// set comm size needed by this Pair if used with fix rigid
comm_forward = 1;
// keep default behavior of history[i][j] = -history[j][i]
nondefault_history_transfer = 0;
}
/* ---------------------------------------------------------------------- */
@ -413,7 +418,7 @@ void PairGranHookeHistory::init_style()
if (history && fix_history == NULL) {
char dnumstr[16];
sprintf(dnumstr,"%d",3);
sprintf(dnumstr,"%d",size_history);
char **fixarg = new char*[4];
fixarg[0] = (char *) "NEIGH_HISTORY";
fixarg[1] = (char *) "all";

2
src/GRANULAR/pair_gran_hooke_history.h
View File

@ -54,6 +54,8 @@ class PairGranHookeHistory : public Pair {
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;
int size_history;
class FixNeighHistory *fix_history;
// storage of rigid body masses for use in granular interactions

1741
src/GRANULAR/pair_granular.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

114
src/GRANULAR/pair_granular.h Normal file
View File

@ -0,0 +1,114 @@
/* ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(granular,PairGranular)
#else
#ifndef LMP_PAIR_GRANULAR_H
#define LMP_PAIR_GRANULAR_H
#include "pair.h"
namespace LAMMPS_NS {
class PairGranular : public Pair {
public:
PairGranular(class LAMMPS *);
~PairGranular();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void reset_dt();
double single(int, int, int, int, double, double, double, double &);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double memory_usage();
protected:
double dt;
int freeze_group_bit;
int use_history;
int neighprev;
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;
double **cut;
class FixNeighHistory *fix_history;
// storage of rigid body masses for use in granular interactions
class Fix *fix_rigid; // ptr to rigid body fix, NULL if none
double *mass_rigid; // rigid mass for owned+ghost atoms
int nmax; // allocated size of mass_rigid
void allocate();
void transfer_history(double*, double*);
private:
int size_history;
int *history_transfer_factors;
// model choices
int **normal_model, **damping_model;
int **tangential_model, **roll_model, **twist_model;
// history flags
int normal_history, tangential_history, roll_history, twist_history;
// indices of history entries
int normal_history_index;
int tangential_history_index;
int roll_history_index;
int twist_history_index;
// per-type material coefficients
double **Emod, **poiss, **Gmod;
// per-type coefficients, set in pair coeff command
double ***normal_coeffs;
double ***tangential_coeffs;
double ***roll_coeffs;
double ***twist_coeffs;
// optional user-specified global cutoff, per-type user-specified cutoffs
double **cutoff_type;
double cutoff_global;
double mix_stiffnessE(double, double, double, double);
double mix_stiffnessG(double, double, double, double);
double mix_geom(double, double);
double pulloff_distance(double, double, int, int);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

2
src/MAKE/Makefile.mpi
View File

@ -12,7 +12,7 @@ SHFLAGS = -fPIC
DEPFLAGS = -M
LINK = mpicxx
LINKFLAGS = -g -O
LINKFLAGS = -g -O3
LIB =
SIZE = size

3
src/Purge.list
View File

@ -25,6 +25,9 @@ style_ntopo.h
# other auto-generated files
lmpinstalledpkgs.h
lmpgitversion.h
# removed on 15 March 2019
fix_wall_gran_omp.h
fix_wall_gran_omp.cpp
# renamed on 7 January 2019
pair_lebedeva.cpp
pair_lebedeva.h

186
src/USER-OMP/fix_wall_gran_omp.cpp
View File

@ -1,186 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include <cmath>
#include "fix_wall_gran_omp.h"
#include "atom.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY};
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
FixWallGranOMP::FixWallGranOMP(LAMMPS *lmp, int narg, char **arg) :
FixWallGran(lmp, narg, arg) { }
/* ---------------------------------------------------------------------- */
void FixWallGranOMP::post_force(int /* vflag */)
{
double vwall[3];
// if just reneighbored:
// update rigid body masses for owned atoms if using FixRigid
// body[i] = which body atom I is in, -1 if none
// mass_body = mass of each rigid body
if (neighbor->ago == 0 && fix_rigid) {
int tmp;
const int * const body = (const int *) fix_rigid->extract("body",tmp);
double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
if (atom->nmax > nmax) {
memory->destroy(mass_rigid);
nmax = atom->nmax;
memory->create(mass_rigid,nmax,"wall/gran:mass_rigid");
}
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
else mass_rigid[i] = 0.0;
}
}
// set position of wall to initial settings and velocity to 0.0
// if wiggle or shear, set wall position and velocity accordingly
double wlo = lo;
double whi = hi;
vwall[0] = vwall[1] = vwall[2] = 0.0;
if (wiggle) {
double arg = omega * (update->ntimestep - time_origin) * dt;
if (wallstyle == axis) {
wlo = lo + amplitude - amplitude*cos(arg);
whi = hi + amplitude - amplitude*cos(arg);
}
vwall[axis] = amplitude*omega*sin(arg);
} else if (wshear) vwall[axis] = vshear;
// loop over all my atoms
// rsq = distance from wall
// dx,dy,dz = signed distance from wall
// for rotating cylinder, reset vwall based on particle position
// skip atom if not close enough to wall
// if wall was set to NULL, it's skipped since lo/hi are infinity
// compute force and torque on atom if close enough to wall
// via wall potential matched to pair potential
// set shear if pair potential stores history
double * const * const x = atom->x;
double * const * const v = atom->v;
double * const * const f = atom->f;
double * const * const omega = atom->omega;
double * const * const torque = atom->torque;
double * const radius = atom->radius;
double * const rmass = atom->rmass;
const int * const mask = atom->mask;
const int nlocal = atom->nlocal;
shearupdate = (update->setupflag) ? 0 : 1;
int i;
#if defined(_OPENMP)
#pragma omp parallel for private(i) default(none) firstprivate(vwall,wlo,whi)
#endif
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double dx,dy,dz,del1,del2,delxy,delr,rsq;
double rwall = 0.0;
dx = dy = dz = 0.0;
if (wallstyle == XPLANE) {
del1 = x[i][0] - wlo;
del2 = whi - x[i][0];
if (del1 < del2) dx = del1;
else dx = -del2;
} else if (wallstyle == YPLANE) {
del1 = x[i][1] - wlo;
del2 = whi - x[i][1];
if (del1 < del2) dy = del1;
else dy = -del2;
} else if (wallstyle == ZPLANE) {
del1 = x[i][2] - wlo;
del2 = whi - x[i][2];
if (del1 < del2) dz = del1;
else dz = -del2;
} else if (wallstyle == ZCYLINDER) {
delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);
delr = cylradius - delxy;
if (delr > radius[i]) {
dz = cylradius;
rwall = 0.0;
} else {
dx = -delr/delxy * x[i][0];
dy = -delr/delxy * x[i][1];
rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius;
if (wshear && axis != 2) {
vwall[0] += vshear * x[i][1]/delxy;
vwall[1] += -vshear * x[i][0]/delxy;
vwall[2] = 0.0;
}
}
}
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) {
if (history)
for (int j = 0; j < sheardim; j++)
shearone[i][j] = 0.0;
} else {
// meff = effective mass of sphere
// if I is part of rigid body, use body mass
double meff = rmass[i];
if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i];
// invoke sphere/wall interaction
if (pairstyle == HOOKE)
hooke(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff);
else if (pairstyle == HOOKE_HISTORY)
hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
else if (pairstyle == HERTZ_HISTORY)
hertz_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
else if (pairstyle == BONDED_HISTORY)
bonded_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixWallGranOMP::post_force_respa(int vflag, int ilevel, int /* iloop */)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}

38
src/USER-OMP/fix_wall_gran_omp.h
View File

@ -1,38 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(wall/gran/omp,FixWallGranOMP)
#else
#ifndef LMP_FIX_WALL_GRAN_OMP_H
#define LMP_FIX_WALL_GRAN_OMP_H
#include "fix_wall_gran.h"
namespace LAMMPS_NS {
class FixWallGranOMP : public FixWallGran {
public:
FixWallGranOMP(class LAMMPS *, int, char **);
virtual void post_force(int);
virtual void post_force_respa(int, int, int);
};
}
#endif
#endif

10
src/fix_neigh_history.cpp
View File

@ -408,7 +408,9 @@ void FixNeighHistory::pre_exchange_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues,jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}
}
@ -520,7 +522,9 @@ void FixNeighHistory::pre_exchange_no_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues, jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}
}
@ -604,7 +608,7 @@ void FixNeighHistory::post_neighbor()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
rflag = sbmask(j);
rflag = sbmask(j) | pair->beyond_contact;
j &= NEIGHMASK;
jlist[jj] = j;

3
src/pair.cpp
View File

@ -103,6 +103,9 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
num_tally_compute = 0;
list_tally_compute = NULL;
nondefault_history_transfer = 0;
beyond_contact = 0;
// KOKKOS per-fix data masks
execution_space = Host;

4
src/pair.h
View File

@ -98,6 +98,8 @@ class Pair : protected Pointers {
enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options
int beyond_contact, nondefault_history_transfer; // for granular styles
// KOKKOS host/device flag and data masks
ExecutionSpace execution_space;
@ -180,6 +182,7 @@ class Pair : protected Pointers {
virtual void min_xf_pointers(int, double **, double **) {}
virtual void min_xf_get(int) {}
virtual void min_x_set(int) {}
virtual void transfer_history(double *, double*) {}
// management of callbacks to be run from ev_tally()
@ -202,6 +205,7 @@ class Pair : protected Pointers {
double tabinner; // inner cutoff for Coulomb table
double tabinner_disp; // inner cutoff for dispersion table
public:
// custom data type for accessing Coulomb tables