lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

cosmetic reformatting of new GRANULAR files

This commit is contained in:
Steve Plimpton 2019-03-27 16:58:14 -06:00
parent 2bac365081
commit 5210c4c3a4
12 changed files with 1086 additions and 825 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
doc/src/fix_wall_gran.txt
View File

@ -59,31 +59,32 @@ 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 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.
"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 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
"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.
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 wall/particle coefficients than for particle/particle
@ -121,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:

59
doc/src/fix_wall_gran_region.txt
View File

@ -48,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)
@ -129,15 +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 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.
"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
@ -150,17 +151,17 @@ 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 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
"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.
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
@ -169,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
@ -185,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

571
doc/src/pair_granular.txt
View File

@ -19,8 +19,7 @@ pair_style granular command :h3
pair_style granular cutoff :pre
cutoff = global cutoff (optional). See discussion below. :l
:ule
cutoff = global cutoff (optional). See discussion below. :l
[Examples:]
@ -44,29 +43,39 @@ pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sd
[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.
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.
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:
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\))
@ -74,137 +83,165 @@ for normal contact models and their required arguments are:
{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.
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:
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}.
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}.
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\}\)
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:
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:
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:
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.
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:
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\}\).
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:
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.
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:
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},
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:
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}).
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:
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:
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 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:
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\}
@ -212,24 +249,24 @@ The total normal force is computed as the sum of the elastic and damping compone
: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:
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
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:
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\}
@ -241,41 +278,52 @@ The tangential damping force \(\mathbf\{F\}_\mathrm\{t,damp\}\) is given by:
\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):
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 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:
({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:
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}.
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
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:
@ -284,49 +332,55 @@ For {tangential linear_history}, the tangential force is given by:
\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:
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:
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.
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):
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:
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
@ -343,72 +397,90 @@ option by an additional factor of {a}, the radius of the contact region. The tan
\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:
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\}\}\):
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.
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:
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:
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:
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:
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.
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:
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\}
@ -420,9 +492,10 @@ but acts only to induce an equal and opposite torque on each particle, according
: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:
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\}\)
@ -430,36 +503,42 @@ centers. The options currently supported are:
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:
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:
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:
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 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):
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
@ -485,19 +564,25 @@ Finally, the twisting torque on each particle is given by:
: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.
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
@ -529,20 +614,21 @@ 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:
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:
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\}
@ -610,57 +696,70 @@ compute depend on atom velocities. See the
[Default:]
For the {pair_coeff} settings: {damping viscoelastic}, {rolling none}, {twisting none}
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.
[(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(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.
[(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.
[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.
[(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.
[(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.
[(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.
[(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.
[(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
[(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.
[(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.
[(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

746
src/GRANULAR/fix_wall_gran.cpp
View File

@ -124,21 +124,27 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
roll_model = twist_model = NONE;
while (iarg < narg) {
if (strcmp(arg[iarg], "hooke") == 0) {
if (iarg + 2 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hooke option");
if (iarg + 2 >= narg)
error->all(FLERR,"Illegal fix wall/gran command, "
"not enough parameters provided for Hooke option");
normal_model = NORMAL_HOOKE;
normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping
iarg += 3;
} else if (strcmp(arg[iarg], "hertz") == 0) {
int num_coeffs = 2;
if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option");
if (iarg + num_coeffs >= narg)
error->all(FLERR,"Illegal fix wall/gran command, "
"not enough parameters provided for Hertz option");
normal_model = NORMAL_HERTZ;
normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping
iarg += num_coeffs+1;
} else if (strcmp(arg[iarg], "hertz/material") == 0) {
int num_coeffs = 3;
if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option");
if (iarg + num_coeffs >= narg)
error->all(FLERR,"Illegal fix wall/gran command, "
"not enough parameters provided for Hertz option");
normal_model = HERTZ_MATERIAL;
Emod = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping
@ -147,7 +153,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
normal_coeffs[2] = poiss;
iarg += num_coeffs+1;
} else if (strcmp(arg[iarg], "dmt") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal fix wall/gran command, "
"not enough parameters provided for Hertz option");
normal_model = DMT;
Emod = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping
@ -157,7 +165,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
iarg += 5;
} else if (strcmp(arg[iarg], "jkr") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for JKR option");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal wall/gran command, "
"not enough parameters provided for JKR option");
beyond_contact = 1;
normal_model = JKR;
Emod = force->numeric(FLERR,arg[iarg+1]); //E
@ -168,7 +178,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
iarg += 5;
} else if (strcmp(arg[iarg], "damping") == 0) {
if (iarg+1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters provided for damping model");
if (iarg+1 >= narg)
error->all(FLERR, "Illegal wall/gran command, "
"not enough parameters provided for damping model");
if (strcmp(arg[iarg+1], "velocity") == 0) {
damping_model = VELOCITY;
iarg += 1;
@ -178,58 +190,80 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg+1], "tsuji") == 0) {
damping_model = TSUJI;
iarg += 1;
} else error->all(FLERR, "Illegal wall/gran command, unrecognized damping model");
} else error->all(FLERR, "Illegal wall/gran command, "
"unrecognized damping model");
iarg += 1;
} else if (strcmp(arg[iarg], "tangential") == 0) {
if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword");
if (iarg + 1 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"must specify tangential model after tangential keyword");
if (strcmp(arg[iarg+1], "linear_nohistory") == 0) {
if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (iarg + 3 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model");
tangential_model = TANGENTIAL_NOHISTORY;
tangential_coeffs[0] = 0;
tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //gammat
tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff.
// gammat and friction coeff
tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+2]);
tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
(strcmp(arg[iarg+1], "mindlin") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (strcmp(arg[iarg+1], "linear_history") == 0) tangential_model = TANGENTIAL_HISTORY;
else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model = TANGENTIAL_MINDLIN_RESCALE;
if ((tangential_model == TANGENTIAL_MINDLIN || tangential_model == TANGENTIAL_MINDLIN_RESCALE) &&
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model");
if (strcmp(arg[iarg+1], "linear_history") == 0)
tangential_model = TANGENTIAL_HISTORY;
else if (strcmp(arg[iarg+1], "mindlin") == 0)
tangential_model = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
tangential_model = TANGENTIAL_MINDLIN_RESCALE;
if ((tangential_model == TANGENTIAL_MINDLIN ||
tangential_model == TANGENTIAL_MINDLIN_RESCALE) &&
(strcmp(arg[iarg+2], "NULL") == 0)) {
if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) {
error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties");
error->all(FLERR, "NULL setting for Mindlin tangential "
"stiffness requires a normal contact model "
"that specifies material properties");
}
tangential_coeffs[0] = 4*(2-poiss)*(1+poiss)/Emod;
} else {
tangential_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt
}
tangential_history = 1;
tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
// gammat and friction coeff
tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]);
tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else {
error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized");
error->all(FLERR, "Illegal pair_coeff command, "
"tangential model not recognized");
}
} else if (strcmp(arg[iarg], "rolling") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters");
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal wall/gran command, not enough parameters");
if (strcmp(arg[iarg+1], "none") == 0) {
roll_model = ROLL_NONE;
iarg += 2;
} else if (strcmp(arg[iarg+1], "sds") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for rolling model");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal wall/gran command, "
"not enough parameters provided for rolling model");
roll_model = ROLL_SDS;
roll_history = 1;
roll_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kR
roll_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR
roll_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff.
// kR, gammaR, rolling friction coeff
roll_coeffs[0] = force->numeric(FLERR,arg[iarg+2]);
roll_coeffs[1] = force->numeric(FLERR,arg[iarg+3]);
roll_coeffs[2] = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else {
error->all(FLERR, "Illegal wall/gran command, rolling friction model not recognized");
error->all(FLERR, "Illegal wall/gran command, "
"rolling friction model not recognized");
}
} else if (strcmp(arg[iarg], "twisting") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters");
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal wall/gran command, not enough parameters");
if (strcmp(arg[iarg+1], "none") == 0) {
twist_model = TWIST_NONE;
iarg += 2;
@ -238,7 +272,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
twist_history = 1;
iarg += 2;
} else if (strcmp(arg[iarg+1], "sds") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for twist model");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal wall/gran command, "
"not enough parameters provided for twist model");
twist_model = TWIST_SDS;
twist_history = 1;
twist_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt
@ -246,7 +282,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
twist_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
} else {
error->all(FLERR, "Illegal wall/gran command, twisting friction model not recognized");
error->all(FLERR, "Illegal wall/gran command, "
"twisting friction model not recognized");
}
} else if (strcmp(arg[iarg], "xplane") == 0 ||
strcmp(arg[iarg], "yplane") == 0 ||
@ -472,8 +509,6 @@ void FixWallGran::init()
27.467*pow(cor,4)-18.022*pow(cor,5)+
4.8218*pow(cor,6);
}
}
/* ---------------------------------------------------------------------- */
@ -616,7 +651,8 @@ void FixWallGran::post_force(int /*vflag*/)
else {
if (pairstyle == GRANULAR && normal_model == JKR && use_history) {
if ((history_one[i][0] == 0) && (rsq > radius[i]*radius[i])) {
// Particles have not contacted yet, and are outside of contact distance
// Particles have not contacted yet,
// and are outside of contact distance
for (j = 0; j < size_history; j++)
history_one[i][j] = 0.0;
continue;
@ -824,7 +860,8 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
history[1] += vtr2*dt;
history[2] += vtr3*dt;
}
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]);
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
// rotate shear displacements
@ -954,7 +991,8 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
history[1] += vtr2*dt;
history[2] += vtr3*dt;
}
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]);
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
// rotate history displacements
@ -1015,372 +1053,382 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
torque[2] -= radius*tor3;
}
/* ---------------------------------------------------------------------- */
void FixWallGran::granular(double rsq, double dx, double dy, double dz,
double *vwall, double rwall, double *v,
double *f, double *omega, double *torque,
double radius, double meff, double *history,
double *contact)
double *vwall, double rwall, double *v,
double *f, double *omega, double *torque,
double radius, double meff, double *history,
double *contact)
{
double fx,fy,fz,nx,ny,nz;
double radsum,r,rinv;
double Reff, delta, dR, dR2;
double fx,fy,fz,nx,ny,nz;
double radsum,r,rinv;
double Reff, delta, dR, dR2;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double knfac, damp_normal, damp_normal_prefactor;
double k_tangential, damp_tangential;
double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;
double fs, fs1, fs2, fs3;
double knfac, damp_normal, damp_normal_prefactor;
double k_tangential, damp_tangential;
double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;
double fs, fs1, fs2, fs3;
double tor1,tor2,tor3;
double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
double tor1,tor2,tor3;
double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
//For JKR
double R2, coh, F_pulloff, a, a2, E;
double t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3;
// for JKR
double R2, coh, F_pulloff, a, a2, E;
double t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3;
//Rolling
double k_roll, damp_roll;
double torroll1, torroll2, torroll3;
double rollmag, rolldotn, scalefac;
double fr, fr1, fr2, fr3;
// rolling
double k_roll, damp_roll;
double torroll1, torroll2, torroll3;
double rollmag, rolldotn, scalefac;
double fr, fr1, fr2, fr3;
//Twisting
double k_twist, damp_twist, mu_twist;
double signtwist, magtwist, magtortwist, Mtcrit;
double tortwist1, tortwist2, tortwist3;
// twisting
double k_twist, damp_twist, mu_twist;
double signtwist, magtwist, magtortwist, Mtcrit;
double tortwist1, tortwist2, tortwist3;
double shrmag,rsht;
double shrmag,rsht;
r = sqrt(rsq);
radsum = rwall + radius;
r = sqrt(rsq);
radsum = rwall + radius;
E = normal_coeffs[0];
E = normal_coeffs[0];
radsum = radius + rwall;
if (rwall == 0) Reff = radius;
else Reff = radius*rwall/(radius+rwall);
radsum = radius + rwall;
if (rwall == 0) Reff = radius;
else Reff = radius*rwall/(radius+rwall);
rinv = 1.0/r;
rinv = 1.0/r;
nx = dx*rinv;
ny = dy*rinv;
nz = dz*rinv;
nx = dx*rinv;
ny = dy*rinv;
nz = dz*rinv;
// relative translational velocity
// relative translational velocity
vr1 = v[0] - vwall[0];
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
vr1 = v[0] - vwall[0];
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
// normal component
// normal component
vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n
vn1 = nx*vnnr;
vn2 = ny*vnnr;
vn3 = nz*vnnr;
vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n
vn1 = nx*vnnr;
vn2 = ny*vnnr;
vn3 = nz*vnnr;
delta = radsum - r;
dR = delta*Reff;
if (normal_model == JKR) {
history[0] = 1.0;
E *= THREEQUARTERS;
R2=Reff*Reff;
coh = normal_coeffs[3];
dR2 = dR*dR;
t0 = coh*coh*R2*R2*E;
t1 = PI27SQ*t0;
t2 = 8*dR*dR2*E*E*E;
t3 = 4*dR2*E;
sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues
t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1));
t5 = t3/t4 + t4/E;
sqrt2 = MAX(0, 2*dR + t5);
t6 = sqrt(sqrt2);
sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6));
a = INVROOT6*(t6 + sqrt(sqrt3));
a2 = a*a;
knfac = normal_coeffs[0]*a;
Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
}
else{
knfac = E; //Hooke
a = sqrt(dR);
if (normal_model != HOOKE) {
Fne *= a;
knfac *= a;
}
Fne = knfac*delta;
if (normal_model == DMT)
Fne -= 4*MY_PI*normal_coeffs[3]*Reff;
}
delta = radsum - r;
dR = delta*Reff;
if (normal_model == JKR) {
history[0] = 1.0;
E *= THREEQUARTERS;
R2=Reff*Reff;
coh = normal_coeffs[3];
dR2 = dR*dR;
t0 = coh*coh*R2*R2*E;
t1 = PI27SQ*t0;
t2 = 8*dR*dR2*E*E*E;
t3 = 4*dR2*E;
sqrt1 = MAX(0, t0*(t1+2*t2)); // in case sqrt(0) < 0 due to precision issues
t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1));
t5 = t3/t4 + t4/E;
sqrt2 = MAX(0, 2*dR + t5);
t6 = sqrt(sqrt2);
sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6));
a = INVROOT6*(t6 + sqrt(sqrt3));
a2 = a*a;
knfac = normal_coeffs[0]*a;
Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
}
else{
knfac = E; //Hooke
a = sqrt(dR);
if (normal_model != HOOKE) {
Fne *= a;
knfac *= a;
}
Fne = knfac*delta;
if (normal_model == DMT)
Fne -= 4*MY_PI*normal_coeffs[3]*Reff;
}
if (damping_model == VELOCITY) {
damp_normal = 1;
}
else if (damping_model == VISCOELASTIC) {
damp_normal = a*meff;
}
else if (damping_model == TSUJI) {
damp_normal = sqrt(meff*knfac);
}
if (damping_model == VELOCITY) {
damp_normal = 1;
}
else if (damping_model == VISCOELASTIC) {
damp_normal = a*meff;
}
else if (damping_model == TSUJI) {
damp_normal = sqrt(meff*knfac);
}
damp_normal_prefactor = normal_coeffs[1]*damp_normal;
Fdamp = -damp_normal_prefactor*vnnr;
damp_normal_prefactor = normal_coeffs[1]*damp_normal;
Fdamp = -damp_normal_prefactor*vnnr;
Fntot = Fne + Fdamp;
Fntot = Fne + Fdamp;
//****************************************
//Tangential force, including history effects
//****************************************
//****************************************
// tangential force, including history effects
//****************************************
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
// relative rotational velocity
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
// relative tangential velocities
vtr1 = vt1 - (nz*wr2-ny*wr3);
vtr2 = vt2 - (nx*wr3-nz*wr1);
vtr3 = vt3 - (ny*wr1-nx*wr2);
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
// relative tangential velocities
vtr1 = vt1 - (nz*wr2-ny*wr3);
vtr2 = vt2 - (nx*wr3-nz*wr1);
vtr3 = vt3 - (ny*wr1-nx*wr2);
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
if (normal_model == JKR) {
F_pulloff = 3*M_PI*coh*Reff;
Fncrit = fabs(Fne + 2*F_pulloff);
}
else if (normal_model == DMT) {
F_pulloff = 4*M_PI*coh*Reff;
Fncrit = fabs(Fne + 2*F_pulloff);
}
else{
Fncrit = fabs(Fntot);
}
if (normal_model == JKR) {
F_pulloff = 3*M_PI*coh*Reff;
Fncrit = fabs(Fne + 2*F_pulloff);
}
else if (normal_model == DMT) {
F_pulloff = 4*M_PI*coh*Reff;
Fncrit = fabs(Fne + 2*F_pulloff);
}
else{
Fncrit = fabs(Fntot);
}
//------------------------------
//Tangential forces
//------------------------------
k_tangential = tangential_coeffs[0];
damp_tangential = tangential_coeffs[1]*damp_normal_prefactor;
//------------------------------
// tangential forces
//------------------------------
int thist0 = tangential_history_index;
int thist1 = thist0 + 1;
int thist2 = thist1 + 1;
k_tangential = tangential_coeffs[0];
damp_tangential = tangential_coeffs[1]*damp_normal_prefactor;
if (tangential_history) {
if (tangential_model == TANGENTIAL_MINDLIN) {
k_tangential *= a;
}
else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a;
if (a < history[3]) { //On unloading, rescale the shear displacements
double factor = a/history[thist2+1];
history[thist0] *= factor;
history[thist1] *= factor;
history[thist2] *= factor;
}
}
shrmag = sqrt(history[thist0]*history[thist0] + history[thist1]*history[thist1] +
history[thist2]*history[thist2]);
int thist0 = tangential_history_index;
int thist1 = thist0 + 1;
int thist2 = thist1 + 1;
// Rotate and update displacements.
// See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (history_update) {
rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz;
if (fabs(rsht) < EPSILON) rsht = 0;
if (rsht > 0) {
scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash!
history[thist0] -= rsht*nx;
history[thist1] -= rsht*ny;
history[thist2] -= rsht*nz;
//Also rescale to preserve magnitude
history[thist0] *= scalefac;
history[thist1] *= scalefac;
history[thist2] *= scalefac;
}
//Update history
history[thist0] += vtr1*dt;
history[thist1] += vtr2*dt;
history[thist2] += vtr3*dt;
}
if (tangential_history) {
if (tangential_model == TANGENTIAL_MINDLIN) {
k_tangential *= a;
}
else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a;
if (a < history[3]) { //On unloading, rescale the shear displacements
double factor = a/history[thist2+1];
history[thist0] *= factor;
history[thist1] *= factor;
history[thist2] *= factor;
}
}
shrmag = sqrt(history[thist0]*history[thist0] +
history[thist1]*history[thist1] +
history[thist2]*history[thist2]);
// tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3;
// rotate and update displacements.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (history_update) {
rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz;
if (fabs(rsht) < EPSILON) rsht = 0;
if (rsht > 0) {
// if rhst == shrmag, contacting pair has rotated 90 deg in one step,
// in which case you deserve a crash!
scalefac = shrmag/(shrmag - rsht);
history[thist0] -= rsht*nx;
history[thist1] -= rsht*ny;
history[thist2] -= rsht*nz;
// also rescale to preserve magnitude
history[thist0] *= scalefac;
history[thist1] *= scalefac;
history[thist2] *= scalefac;
}
// update history
history[thist0] += vtr1*dt;
history[thist1] += vtr2*dt;
history[thist2] += vtr3*dt;
}
// rescale frictional displacements and forces if needed
Fscrit = tangential_coeffs[2] * Fncrit;
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
if (fs > Fscrit) {
if (shrmag != 0.0) {
history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1);
history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2);
history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3);
fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs;
} else fs1 = fs2 = fs3 = 0.0;
}
}
else{ //Classic pair gran/hooke (no history)
fs = meff*damp_tangential*vrel;
if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
else Ft = 0.0;
fs1 = -Ft*vtr1;
fs2 = -Ft*vtr2;
fs3 = -Ft*vtr3;
}
// tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3;
//****************************************
// Rolling resistance
//****************************************
// rescale frictional displacements and forces if needed
Fscrit = tangential_coeffs[2] * Fncrit;
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
if (fs > Fscrit) {
if (shrmag != 0.0) {
history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs +
damp_tangential*vtr1);
history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs +
damp_tangential*vtr2);
history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs +
damp_tangential*vtr3);
fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs;
} else fs1 = fs2 = fs3 = 0.0;
}
} else { // classic pair gran/hooke (no history)
fs = meff*damp_tangential*vrel;
if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
else Ft = 0.0;
fs1 = -Ft*vtr1;
fs2 = -Ft*vtr2;
fs3 = -Ft*vtr3;
}
if (roll_model != ROLL_NONE) {
relrot1 = omega[0];
relrot2 = omega[1];
relrot3 = omega[2];
//****************************************
// rolling resistance
//****************************************
// rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// This is different from the Marshall papers, which use the Bagi/Kuhn formulation
// for rolling velocity (see Wang et al for why the latter is wrong)
vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1;
vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2;
vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3;
if (roll_model != ROLL_NONE) {
relrot1 = omega[0];
relrot2 = omega[1];
relrot3 = omega[2];
int rhist0 = roll_history_index;
int rhist1 = rhist0 + 1;
int rhist2 = rhist1 + 1;
// rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// This is different from the Marshall papers,
// which use the Bagi/Kuhn formulation
// for rolling velocity (see Wang et al for why the latter is wrong)
vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1;
vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2;
vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3;
// Rolling displacement
rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
int rhist0 = roll_history_index;
int rhist1 = rhist0 + 1;
int rhist2 = rhist1 + 1;
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
// rolling displacement
rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
if (history_update) {
if (fabs(rolldotn) < EPSILON) rolldotn = 0;
if (rolldotn > 0) { //Rotate into tangential plane
scalefac = rollmag/(rollmag - rolldotn);
history[rhist0] -= rolldotn*nx;
history[rhist1] -= rolldotn*ny;
history[rhist2] -= rolldotn*nz;
//Also rescale to preserve magnitude
history[rhist0] *= scalefac;
history[rhist1] *= scalefac;
history[rhist2] *= scalefac;
}
history[rhist0] += vrl1*dt;
history[rhist1] += vrl2*dt;
history[rhist2] += vrl3*dt;
}
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
k_roll = roll_coeffs[0];
damp_roll = roll_coeffs[1];
fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
if (history_update) {
if (fabs(rolldotn) < EPSILON) rolldotn = 0;
if (rolldotn > 0) { // rotate into tangential plane
scalefac = rollmag/(rollmag - rolldotn);
history[rhist0] -= rolldotn*nx;
history[rhist1] -= rolldotn*ny;
history[rhist2] -= rolldotn*nz;
// also rescale to preserve magnitude
history[rhist0] *= scalefac;
history[rhist1] *= scalefac;
history[rhist2] *= scalefac;
}
history[rhist0] += vrl1*dt;
history[rhist1] += vrl2*dt;
history[rhist2] += vrl3*dt;
}
// rescale frictional displacements and forces if needed
Frcrit = roll_coeffs[2] * Fncrit;
k_roll = roll_coeffs[0];
damp_roll = roll_coeffs[1];
fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
if (fr > Frcrit) {
if (rollmag != 0.0) {
history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1);
history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2);
history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3);
fr1 *= Frcrit/fr;
fr2 *= Frcrit/fr;
fr3 *= Frcrit/fr;
} else fr1 = fr2 = fr3 = 0.0;
}
}
// rescale frictional displacements and forces if needed
Frcrit = roll_coeffs[2] * Fncrit;
//****************************************
// Twisting torque, including history effects
//****************************************
if (twist_model != TWIST_NONE) {
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
if (twist_model == TWIST_MARSHALL) {
k_twist = 0.5*k_tangential*a*a;; //eq 32 of Marshall paper
damp_twist = 0.5*damp_tangential*a*a;
mu_twist = TWOTHIRDS*a*tangential_coeffs[2];
}
else{
k_twist = twist_coeffs[0];
damp_twist = twist_coeffs[1];
mu_twist = twist_coeffs[2];
}
if (history_update) {
history[twist_history_index] += magtwist*dt;
}
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30)
signtwist = (magtwist > 0) - (magtwist < 0);
Mtcrit = mu_twist*Fncrit;//critical torque (eq 44)
if (fabs(magtortwist) > Mtcrit) {
history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist);
magtortwist = -Mtcrit * signtwist; //eq 34
}
}
// Apply forces & torques
fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
if (fr > Frcrit) {
if (rollmag != 0.0) {
history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1);
history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2);
history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3);
fr1 *= Frcrit/fr;
fr2 *= Frcrit/fr;
fr3 *= Frcrit/fr;
} else fr1 = fr2 = fr3 = 0.0;
}
}
fx = nx*Fntot + fs1;
fy = ny*Fntot + fs2;
fz = nz*Fntot + fs3;
//****************************************
// twisting torque, including history effects
//****************************************
if (peratom_flag) {
contact[1] = fx;
contact[2] = fy;
contact[3] = fz;
}
if (twist_model != TWIST_NONE) {
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
if (twist_model == TWIST_MARSHALL) {
k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper
damp_twist = 0.5*damp_tangential*a*a;
mu_twist = TWOTHIRDS*a*tangential_coeffs[2];
}
else{
k_twist = twist_coeffs[0];
damp_twist = twist_coeffs[1];
mu_twist = twist_coeffs[2];
}
if (history_update) {
history[twist_history_index] += magtwist*dt;
}
// M_t torque (eq 30)
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;
signtwist = (magtwist > 0) - (magtwist < 0);
Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
if (fabs(magtortwist) > Mtcrit) {
history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist -
damp_twist*magtwist);
magtortwist = -Mtcrit * signtwist; // eq 34
}
}
f[0] += fx;
f[1] += fy;
f[2] += fz;
// apply forces & torques
tor1 = ny*fs3 - nz*fs2;
tor2 = nz*fs1 - nx*fs3;
tor3 = nx*fs2 - ny*fs1;
fx = nx*Fntot + fs1;
fy = ny*Fntot + fs2;
fz = nz*Fntot + fs3;
torque[0] -= radius*tor1;
torque[1] -= radius*tor2;
torque[2] -= radius*tor3;
if (peratom_flag) {
contact[1] = fx;
contact[2] = fy;
contact[3] = fz;
}
if (twist_model != TWIST_NONE) {
tortwist1 = magtortwist * nx;
tortwist2 = magtortwist * ny;
tortwist3 = magtortwist * nz;
f[0] += fx;
f[1] += fy;
f[2] += fz;
torque[0] += tortwist1;
torque[1] += tortwist2;
torque[2] += tortwist3;
}
tor1 = ny*fs3 - nz*fs2;
tor2 = nz*fs1 - nx*fs3;
tor3 = nx*fs2 - ny*fs1;
if (roll_model != ROLL_NONE) {
torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr
torroll2 = Reff*(nz*fr1 - nx*fr3);
torroll3 = Reff*(nx*fr2 - ny*fr1);
torque[0] -= radius*tor1;
torque[1] -= radius*tor2;
torque[2] -= radius*tor3;
torque[0] += torroll1;
torque[1] += torroll2;
torque[2] += torroll3;
}
if (twist_model != TWIST_NONE) {
tortwist1 = magtortwist * nx;
tortwist2 = magtortwist * ny;
tortwist3 = magtortwist * nz;
torque[0] += tortwist1;
torque[1] += tortwist2;
torque[2] += tortwist3;
}
if (roll_model != ROLL_NONE) {
torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr
torroll2 = Reff*(nz*fr1 - nx*fr3);
torroll3 = Reff*(nx*fr2 - ny*fr1);
torque[0] += torroll1;
torque[1] += torroll2;
torque[2] += torroll3;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
@ -1389,9 +1437,10 @@ double FixWallGran::memory_usage()
{
int nmax = atom->nmax;
double bytes = 0.0;
if (use_history) bytes += nmax*size_history * sizeof(double); // shear history
if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid
if (peratom_flag) bytes += nmax*size_peratom_cols*sizeof(double); //store contacts
if (use_history) bytes += nmax*size_history * sizeof(double); // shear history
if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid
// store contacts
if (peratom_flag) bytes += nmax*size_peratom_cols*sizeof(double);
return bytes;
}
@ -1401,7 +1450,8 @@ double FixWallGran::memory_usage()
void FixWallGran::grow_arrays(int nmax)
{
if (use_history) memory->grow(history_one,nmax,size_history,"fix_wall_gran:history_one");
if (use_history) memory->grow(history_one,nmax,size_history,
"fix_wall_gran:history_one");
if (peratom_flag) {
memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");
}

17
src/GRANULAR/fix_wall_gran.h
View File

@ -66,26 +66,24 @@ class FixWallGran : public Fix {
bigint time_origin;
double kn,kt,gamman,gammat,xmu;
//For granular
//Model choices
// for granular model choices
int normal_model, damping_model;
int tangential_model, roll_model, twist_model;
int beyond_contact;
//History flags
// history flags
int normal_history, tangential_history, roll_history, twist_history;
//Indices of history entries
// indices of history entries
int normal_history_index;
int tangential_history_index;
int roll_history_index;
int twist_history_index;
//Material coefficients
// material coefficients
double Emod, poiss, Gmod;
//Contact model coefficients
// contact model coefficients
double normal_coeffs[4];
double tangential_coeffs[3];
double roll_coeffs[3];
@ -97,7 +95,7 @@ class FixWallGran : public Fix {
char *idregion;
int use_history; // if particle/wall interaction stores history
int history_update; // flag for whether shear history is updated
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
@ -110,7 +108,8 @@ class FixWallGran : public Fix {
double *mass_rigid; // rigid mass for owned+ghost atoms
int nmax; // allocated size of mass_rigid
// Store particle interactions
// store particle interactions
int store;
};

19
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -47,8 +47,9 @@ enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR};
/* ---------------------------------------------------------------------- */
FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), ncontact(NULL),
walls(NULL), history_many(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;
@ -190,7 +191,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
if (!region->match(x[i][0],x[i][1],x[i][2])) continue;
if (pairstyle == GRANULAR && normal_model == JKR){
nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]+pulloff_distance(radius[i]));
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]);
@ -278,8 +280,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
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);
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);
}
}
}
@ -364,11 +367,11 @@ void FixWallGranRegion::grow_arrays(int nmax)
if (use_history) {
memory->grow(ncontact,nmax,"fix_wall_gran:ncontact");
memory->grow(walls,nmax,tmax,"fix_wall_gran:walls");
memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many");
memory->grow(history_many,nmax,tmax,size_history,
"fix_wall_gran:history_many");
}
if (peratom_flag){
if (peratom_flag)
memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");
}
}
/* ----------------------------------------------------------------------

4
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -59,7 +59,9 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
comm_forward = 1;
nondefault_history_transfer = 0; //keep default behavior of history[i][j] = -history[j][i]
// keep default behavior of history[i][j] = -history[j][i]
nondefault_history_transfer = 0;
}
/* ---------------------------------------------------------------------- */

3
src/GRANULAR/pair_gran_hooke_history.h
View File

@ -26,6 +26,8 @@ namespace LAMMPS_NS {
class PairGranHookeHistory : public Pair {
public:
int nondefault_history_transfer;
PairGranHookeHistory(class LAMMPS *);
virtual ~PairGranHookeHistory();
virtual void compute(int, int);
@ -42,7 +44,6 @@ class PairGranHookeHistory : public Pair {
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double memory_usage();
int nondefault_history_transfer;
protected:
double kn,kt,gamman,gammat,xmu;

419
src/GRANULAR/pair_granular.cpp
View File

@ -11,9 +11,9 @@ See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL)
Leo Silbert (SNL), Gary Grest (SNL)
Contributing authors:
Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL)
Leo Silbert (SNL), Gary Grest (SNL)
----------------------------------------------------------------------- */
#include <cmath>
@ -52,7 +52,8 @@ using namespace MathConst;
enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR};
enum {VELOCITY, VISCOELASTIC, TSUJI};
enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE};
enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE};
enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
enum {ROLL_NONE, ROLL_SDS};
@ -90,10 +91,10 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
nondefault_history_transfer = 0;
tangential_history_index = 0;
roll_history_index = twist_history_index = 0;
}
/* ---------------------------------------------------------------------- */
PairGranular::~PairGranular()
{
delete [] svector;
@ -118,16 +119,17 @@ PairGranular::~PairGranular()
memory->destroy(roll_model);
memory->destroy(twist_model);
delete [] onerad_dynamic;
delete [] onerad_frozen;
delete [] maxrad_dynamic;
delete [] maxrad_frozen;
}
memory->destroy(mass_rigid);
}
/* ---------------------------------------------------------------------- */
void PairGranular::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
@ -147,18 +149,18 @@ void PairGranular::compute(int eflag, int vflag)
double mi,mj,meff;
double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
// For JKR
// for JKR
double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E;
double t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3;
// Rolling
// rolling
double k_roll, damp_roll;
double torroll1, torroll2, torroll3;
double rollmag, rolldotn, scalefac;
double fr, fr1, fr2, fr3;
// Twisting
// twisting
double k_twist, damp_twist, mu_twist;
double signtwist, magtwist, magtortwist, Mtcrit;
double tortwist1, tortwist2, tortwist3;
@ -317,7 +319,8 @@ void PairGranular::compute(int eflag, int vflag)
t1 = PI27SQ*t0;
t2 = 8*dR*dR2*E*E*E;
t3 = 4*dR2*E;
sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues
// in case sqrt(0) < 0 due to precision issues
sqrt1 = MAX(0, t0*(t1+2*t2));
t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1));
t5 = t3/t4 + t4/E;
sqrt2 = MAX(0, 2*dR + t5);
@ -328,7 +331,7 @@ void PairGranular::compute(int eflag, int vflag)
knfac = normal_coeffs[itype][jtype][0]*a;
Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
} else {
knfac = E; //Hooke
knfac = E; // Hooke
Fne = knfac*delta;
a = sqrt(dR);
if (normal_model[itype][jtype] != HOOKE) {
@ -339,7 +342,9 @@ void PairGranular::compute(int eflag, int vflag)
Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
}
//Consider restricting Hooke to only have 'velocity' as an option for damping?
// NOTE: consider restricting Hooke to only have
// 'velocity' as an option for damping?
if (damping_model[itype][jtype] == VELOCITY) {
damp_normal = 1;
} else if (damping_model[itype][jtype] == VISCOELASTIC) {
@ -354,7 +359,7 @@ void PairGranular::compute(int eflag, int vflag)
Fntot = Fne + Fdamp;
//****************************************
//Tangential force, including history effects
// tangential force, including history effects
//****************************************
// tangential component
@ -374,7 +379,7 @@ void PairGranular::compute(int eflag, int vflag)
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
// If any history is needed:
// if any history is needed
if (use_history) {
touch[jj] = 1;
history = &allhistory[size_history*jj];
@ -391,45 +396,51 @@ void PairGranular::compute(int eflag, int vflag)
}
//------------------------------
//Tangential forces
// tangential forces
//------------------------------
k_tangential = tangential_coeffs[itype][jtype][0];
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
damp_tangential = tangential_coeffs[itype][jtype][1] *
damp_normal_prefactor;
if (tangential_history) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) {
k_tangential *= a;
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
} else if (tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a;
if (a < history[3]) { //On unloading, rescale the shear displacements
// on unloading, rescale the shear displacements
if (a < history[3]) {
double factor = a/history[3];
history[0] *= factor;
history[1] *= factor;
history[2] *= factor;
}
}
// Rotate and update displacements.
// See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
// rotate and update displacements.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (historyupdate) {
rsht = history[0]*nx + history[1]*ny + history[2]*nz;
if (fabs(rsht) < EPSILON) rsht = 0;
if (rsht > 0) {
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
scalefac = shrmag/(shrmag - rsht); //if rsht == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash!
// if rsht == shrmag, contacting pair has rotated 90 deg
// in one step, in which case you deserve a crash!
scalefac = shrmag/(shrmag - rsht);
history[0] -= rsht*nx;
history[1] -= rsht*ny;
history[2] -= rsht*nz;
//Also rescale to preserve magnitude
// also rescale to preserve magnitude
history[0] *= scalefac;
history[1] *= scalefac;
history[2] *= scalefac;
}
//Update history
// update history
history[0] += vtr1*dt;
history[1] += vtr2*dt;
history[2] += vtr3*dt;
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) history[3] = a;
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE)
history[3] = a;
}
// tangential forces = history + tangential velocity damping
@ -444,15 +455,18 @@ void PairGranular::compute(int eflag, int vflag)
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
if (shrmag != 0.0) {
history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1);
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2);
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3);
history[0] = -1.0/k_tangential*(Fscrit*fs1/fs +
damp_tangential*vtr1);
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs +
damp_tangential*vtr2);
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs +
damp_tangential*vtr3);
fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs;
} else fs1 = fs2 = fs3 = 0.0;
}
} else{ //Classic pair gran/hooke (no history)
} else { // classic pair gran/hooke (no history)
fs = meff*damp_tangential*vrel;
if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
else Ft = 0.0;
@ -462,7 +476,7 @@ void PairGranular::compute(int eflag, int vflag)
}
//****************************************
// Rolling resistance
// rolling resistance
//****************************************
if (roll_model[itype][jtype] != ROLL_NONE) {
@ -470,12 +484,18 @@ void PairGranular::compute(int eflag, int vflag)
relrot2 = omega[i][1] - omega[j][1];
relrot3 = omega[i][2] - omega[j][2];
// rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// This is different from the Marshall papers, which use the Bagi/Kuhn formulation
// rolling velocity,
// see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// this is different from the Marshall papers,
// which use the Bagi/Kuhn formulation
// for rolling velocity (see Wang et al for why the latter is wrong)
vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1;
vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2;
vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3;
// - 0.5*((radj-radi)/radsum)*vtr1;
// - 0.5*((radj-radi)/radsum)*vtr2;
// - 0.5*((radj-radi)/radsum)*vtr3;
vrl1 = Reff*(relrot2*nz - relrot3*ny);
vrl2 = Reff*(relrot3*nx - relrot1*nz);
vrl3 = Reff*(relrot1*ny - relrot2*nx);
int rhist0 = roll_history_index;
int rhist1 = rhist0 + 1;
@ -484,7 +504,7 @@ void PairGranular::compute(int eflag, int vflag)
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
if (historyupdate) {
if (fabs(rolldotn) < EPSILON) rolldotn = 0;
if (rolldotn > 0) { //Rotate into tangential plane
if (rolldotn > 0) { // rotate into tangential plane
rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
@ -492,7 +512,7 @@ void PairGranular::compute(int eflag, int vflag)
history[rhist0] -= rolldotn*nx;
history[rhist1] -= rolldotn*ny;
history[rhist2] -= rolldotn*nz;
//Also rescale to preserve magnitude
// also rescale to preserve magnitude
history[rhist0] *= scalefac;
history[rhist1] *= scalefac;
history[rhist2] *= scalefac;
@ -528,12 +548,14 @@ void PairGranular::compute(int eflag, int vflag)
}
//****************************************
// Twisting torque, including history effects
// twisting torque, including history effects
//****************************************
if (twist_model[itype][jtype] != TWIST_NONE) {
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
// omega_T (eq 29 of Marshall)
magtwist = relrot1*nx + relrot2*ny + relrot3*nz;
if (twist_model[itype][jtype] == TWIST_MARSHALL) {
k_twist = 0.5*k_tangential*a*a;; //eq 32 of Marshall paper
k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper
damp_twist = 0.5*damp_tangential*a*a;
mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2];
} else {
@ -544,15 +566,18 @@ void PairGranular::compute(int eflag, int vflag)
if (historyupdate) {
history[twist_history_index] += magtwist*dt;
}
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30)
magtortwist = -k_twist*history[twist_history_index] -
damp_twist*magtwist; // M_t torque (eq 30)
signtwist = (magtwist > 0) - (magtwist < 0);
Mtcrit = mu_twist*Fncrit;//critical torque (eq 44)
Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
if (fabs(magtortwist) > Mtcrit) {
history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist);
magtortwist = -Mtcrit * signtwist; //eq 34
history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist -
damp_twist*magtwist);
magtortwist = -Mtcrit * signtwist; // eq 34
}
}
// Apply forces & torques
// apply forces & torques
fx = nx*Fntot + fs1;
fy = ny*Fntot + fs2;
@ -582,7 +607,7 @@ void PairGranular::compute(int eflag, int vflag)
}
if (roll_model[itype][jtype] != ROLL_NONE) {
torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr
torroll1 = Reff*(ny*fr3 - nz*fr2); // n cross fr
torroll2 = Reff*(nz*fr1 - nx*fr3);
torroll3 = Reff*(nx*fr2 - ny*fr1);
@ -619,9 +644,8 @@ void PairGranular::compute(int eflag, int vflag)
}
}
/* ----------------------------------------------------------------------
allocate all arrays
allocate all arrays
------------------------------------------------------------------------- */
void PairGranular::allocate()
@ -657,7 +681,7 @@ void PairGranular::allocate()
}
/* ----------------------------------------------------------------------
global settings
global settings
------------------------------------------------------------------------- */
void PairGranular::settings(int narg, char **arg)
@ -665,7 +689,7 @@ void PairGranular::settings(int narg, char **arg)
if (narg == 1) {
cutoff_global = force->numeric(FLERR,arg[0]);
} else {
cutoff_global = -1; //Will be set based on particle sizes, model choice
cutoff_global = -1; // will be set based on particle sizes, model choice
}
normal_history = tangential_history = 0;
@ -705,45 +729,57 @@ void PairGranular::coeff(int narg, char **arg)
int iarg = 2;
while (iarg < narg) {
if (strcmp(arg[iarg], "hooke") == 0) {
if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option");
if (iarg + 2 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for Hooke option");
normal_model_one = HOOKE;
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // kn
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping
iarg += 3;
} else if (strcmp(arg[iarg], "hertz") == 0) {
int num_coeffs = 2;
if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option");
if (iarg + num_coeffs >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for Hertz option");
normal_model_one = HERTZ;
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // kn
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping
iarg += num_coeffs+1;
} else if (strcmp(arg[iarg], "hertz/material") == 0) {
int num_coeffs = 3;
if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz/material option");
if (iarg + num_coeffs >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for Hertz/material option");
normal_model_one = HERTZ_MATERIAL;
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio
iarg += num_coeffs+1;
} else if (strcmp(arg[iarg], "dmt") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for Hertz option");
normal_model_one = DMT;
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); // cohesion
iarg += 5;
} else if (strcmp(arg[iarg], "jkr") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for JKR option");
beyond_contact = 1;
normal_model_one = JKR;
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E
normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping
normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio
normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); // cohesion
iarg += 5;
} else if (strcmp(arg[iarg], "damping") == 0) {
if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model");
if (iarg+1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, "
"not enough parameters provided for damping model");
if (strcmp(arg[iarg+1], "velocity") == 0) {
damping_model_one = VELOCITY;
iarg += 1;
@ -753,58 +789,80 @@ void PairGranular::coeff(int narg, char **arg)
} else if (strcmp(arg[iarg+1], "tsuji") == 0) {
damping_model_one = TSUJI;
iarg += 1;
} else error->all(FLERR, "Illegal pair_coeff command, unrecognized damping model");
} else error->all(FLERR, "Illegal pair_coeff command, "
"unrecognized damping model");
iarg += 1;
} else if (strcmp(arg[iarg], "tangential") == 0) {
if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword");
if (iarg + 1 >= narg)
error->all(FLERR,"Illegal pair_coeff command, must specify "
"tangential model after tangential keyword");
if (strcmp(arg[iarg+1], "linear_nohistory") == 0) {
if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (iarg + 3 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model");
tangential_model_one = TANGENTIAL_NOHISTORY;
tangential_coeffs_one[0] = 0;
tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //gammat
tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff.
// gammat and friction coeff
tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]);
tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
(strcmp(arg[iarg+1], "mindlin") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (strcmp(arg[iarg+1], "linear_history") == 0) tangential_model_one = TANGENTIAL_HISTORY;
else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model_one = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model");
if (strcmp(arg[iarg+1], "linear_history") == 0)
tangential_model_one = TANGENTIAL_HISTORY;
else if (strcmp(arg[iarg+1], "mindlin") == 0)
tangential_model_one = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
tangential_history = 1;
if ((tangential_model_one == TANGENTIAL_MINDLIN || tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) &&
if ((tangential_model_one == TANGENTIAL_MINDLIN ||
tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) &&
(strcmp(arg[iarg+2], "NULL") == 0)) {
if (normal_model_one == HERTZ || normal_model_one == HOOKE) {
error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties");
error->all(FLERR, "NULL setting for Mindlin tangential "
"stiffness requires a normal contact model that "
"specifies material properties");
}
tangential_coeffs_one[0] = -1;
} else {
tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt
tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); // kt
}
tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
// gammat and friction coeff
tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]);
tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else {
error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized");
error->all(FLERR, "Illegal pair_coeff command, "
"tangential model not recognized");
}
} else if (strcmp(arg[iarg], "rolling") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strcmp(arg[iarg+1], "none") == 0) {
roll_model_one = ROLL_NONE;
iarg += 2;
} else if (strcmp(arg[iarg+1], "sds") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for rolling model");
roll_model_one = ROLL_SDS;
roll_history = 1;
roll_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kR
roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR
roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff.
// kR and gammaR and rolling friction coeff
roll_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]);
roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]);
roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else {
error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized");
error->all(FLERR, "Illegal pair_coeff command, "
"rolling friction model not recognized");
}
} else if (strcmp(arg[iarg], "twisting") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strcmp(arg[iarg+1], "none") == 0) {
twist_model_one = TWIST_NONE;
iarg += 2;
@ -813,24 +871,31 @@ void PairGranular::coeff(int narg, char **arg)
twist_history = 1;
iarg += 2;
} else if (strcmp(arg[iarg+1], "sds") == 0) {
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model");
if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for twist model");
twist_model_one = TWIST_SDS;
twist_history = 1;
twist_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt
twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
// kt and gammat and friction coeff
twist_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]);
twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]);
twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else {
error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized");
error->all(FLERR, "Illegal pair_coeff command, "
"twisting friction model not recognized");
}
} else if (strcmp(arg[iarg], "cutoff") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (iarg + 1 >= narg)
error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
cutoff_one = force->numeric(FLERR,arg[iarg+1]);
} else error->all(FLERR, "Illegal pair coeff command");
}
//It is an error not to specify normal or tangential model
if ((normal_model_one < 0) || (tangential_model_one < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model");
// error not to specify normal or tangential model
if ((normal_model_one < 0) || (tangential_model_one < 0))
error->all(FLERR, "Illegal pair coeff command, "
"must specify normal or tangential contact model");
int count = 0;
double damp;
@ -849,7 +914,9 @@ void PairGranular::coeff(int narg, char **arg)
if (normal_model_one != HERTZ && normal_model_one != HOOKE) {
Emod[i][j] = Emod[j][i] = normal_coeffs_one[0];
poiss[i][j] = poiss[j][i] = normal_coeffs_one[2];
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = FOURTHIRDS*mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
FOURTHIRDS*mix_stiffnessE(Emod[i][j],Emod[i][j],
poiss[i][j],poiss[i][j]);
} else {
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0];
}
@ -860,12 +927,15 @@ void PairGranular::coeff(int narg, char **arg)
tangential_model[i][j] = tangential_model[j][i] = tangential_model_one;
if (tangential_coeffs_one[0] == -1) {
tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = 8*mix_stiffnessG(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);
tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] =
8*mix_stiffnessG(Emod[i][j],Emod[i][j],poiss[i][j],poiss[i][j]);
} else {
tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = tangential_coeffs_one[0];
tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] =
tangential_coeffs_one[0];
}
for (int k = 1; k < 3; k++)
tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_one[k];
tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] =
tangential_coeffs_one[k];
roll_model[i][j] = roll_model[j][i] = roll_model_one;
if (roll_model_one != ROLL_NONE)
@ -877,13 +947,13 @@ void PairGranular::coeff(int narg, char **arg)
for (int k = 0; k < 3; k++)
twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k];
cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
@ -902,17 +972,21 @@ void PairGranular::init_style()
if (comm->ghost_velocity == 0)
error->all(FLERR,"Pair granular requires ghost atoms store velocity");
// Determine whether we need a granular neigh list, how large it needs to be
use_history = normal_history || tangential_history || roll_history || twist_history;
// determine whether we need a granular neigh list, how large it needs to be
use_history = normal_history || tangential_history ||
roll_history || twist_history;
// for JKR, will need fix/neigh/history to keep track of touch arrays
//For JKR, will need fix/neigh/history to keep track of touch arrays
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (normal_model[i][j] == JKR) use_history = 1;
size_history = 3*tangential_history + 3*roll_history + twist_history;
//Determine location of tangential/roll/twist histories in array
// determine location of tangential/roll/twist histories in array
if (roll_history) {
if (tangential_history) roll_history_index = 3;
else roll_history_index = 0;
@ -1036,7 +1110,7 @@ void PairGranular::init_style()
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairGranular::init_one(int i, int j)
@ -1051,46 +1125,58 @@ double PairGranular::init_one(int i, int j)
(twist_model[i][i] != twist_model[j][j])) {
char str[512];
sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j);
sprintf(str,"Granular pair style functional forms are different, "
"cannot mix coefficients for types %d and %d. \n"
"This combination must be set explicitly "
"via pair_coeff command",i,j);
error->one(FLERR,str);
}
if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE)
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]);
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]);
else
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]);
normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]);
normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]);
normal_coeffs[i][j][1] = normal_coeffs[j][i][1] =
mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]);
if ((normal_model[i][j] == JKR) || (normal_model[i][j] == DMT))
normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]);
normal_coeffs[i][j][3] = normal_coeffs[j][i][3] =
mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]);
for (int k = 0; k < 3; k++)
tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]);
tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] =
mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]);
if (roll_model[i][j] != ROLL_NONE) {
for (int k = 0; k < 3; k++)
roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]);
roll_coeffs[i][j][k] = roll_coeffs[j][i][k] =
mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]);
}
if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL) {
for (int k = 0; k < 3; k++)
twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]);
twist_coeffs[i][j][k] = twist_coeffs[j][i][k] =
mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]);
}
}
// It is possible that cut[i][j] at this point is still 0.0. This can happen when
// It is possible that cut[i][j] at this point is still 0.0.
// This can happen when
// there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates
// problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size
// To avoid this issue, for cases involving cut[i][j] = 0.0 (possible only
// if there is no current information about radius/cutoff of type i and j).
// we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0.
double pulloff;
if (cutoff_type[i][j] < 0 && cutoff_global < 0) {
if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) ||
((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) ||
((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist
// radius info about both i and j exist
((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) {
cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
pulloff = 0.0;
if (normal_model[i][j] == JKR) {
@ -1105,7 +1191,13 @@ double PairGranular::init_one(int i, int j)
if (normal_model[i][j] == JKR)
pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j);
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff);
} else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list)
} else {
// radius info about either i or j does not exist
// (i.e. not present and not about to get poured;
// set to largest value to not interfere with neighbor list)
double cutmax = 0.0;
for (int k = 1; k <= atom->ntypes; k++) {
cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]);
@ -1123,8 +1215,8 @@ double PairGranular::init_one(int i, int j)
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairGranular::write_restart(FILE *fp)
{
@ -1149,8 +1241,8 @@ void PairGranular::write_restart(FILE *fp)
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairGranular::read_restart(FILE *fp)
{
@ -1189,7 +1281,6 @@ void PairGranular::read_restart(FILE *fp)
}
}
/* ---------------------------------------------------------------------- */
void PairGranular::reset_dt()
@ -1216,18 +1307,17 @@ double PairGranular::single(int i, int j, int itype, int jtype,
double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;
double fs, fs1, fs2, fs3;
//For JKR
// for JKR
double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E;
double delta, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3;
//Rolling
// rolling
double k_roll, damp_roll;
double rollmag;
double fr, fr1, fr2, fr3;
//Twisting
// twisting
double k_twist, damp_twist, mu_twist;
double signtwist, magtwist, magtortwist, Mtcrit;
@ -1338,7 +1428,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
t1 = PI27SQ*t0;
t2 = 8*dR*dR2*E*E*E;
t3 = 4*dR2*E;
sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues
// in case sqrt(0) < 0 due to precision issues
sqrt1 = MAX(0, t0*(t1+2*t2));
t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1));
t5 = t3/t4 + t4/E;
sqrt2 = MAX(0, 2*dR + t5);
@ -1387,7 +1478,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
}
//****************************************
//Tangential force, including history effects
// tangential force, including history effects
//****************************************
// tangential component
@ -1418,7 +1509,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
}
//------------------------------
//Tangential forces
// tangential forces
//------------------------------
k_tangential = tangential_coeffs[itype][jtype][0];
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
@ -1428,7 +1519,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
k_tangential *= a;
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a;
if (a < history[3]) { //On unloading, rescale the shear displacements
// on unloading, rescale the shear displacements
if (a < history[3]) {
double factor = a/history[3];
history[0] *= factor;
history[1] *= factor;
@ -1457,7 +1549,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
fs3 *= Fscrit/fs;
} else fs1 = fs2 = fs3 = 0.0;
}
} else { //Classic pair gran/hooke (no history)
// classic pair gran/hooke (no history)
} else {
fs = meff*damp_tangential*vrel;
if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
else Ft = 0.0;
@ -1467,7 +1561,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
}
//****************************************
// Rolling resistance
// rolling resistance
//****************************************
if (roll_model[itype][jtype] != ROLL_NONE) {
@ -1476,8 +1570,10 @@ double PairGranular::single(int i, int j, int itype, int jtype,
relrot3 = omega[i][2] - omega[j][2];
// rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// This is different from the Marshall papers, which use the Bagi/Kuhn formulation
// this is different from the Marshall papers,
// which use the Bagi/Kuhn formulation
// for rolling velocity (see Wang et al for why the latter is wrong)
vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1;
vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2;
vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3;
@ -1486,7 +1582,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
int rhist1 = rhist0 + 1;
int rhist2 = rhist1 + 1;
// Rolling displacement
// rolling displacement
rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
@ -1512,10 +1608,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
}
//****************************************
// Twisting torque, including history effects
// twisting torque, including history effects
//****************************************
if (twist_model[itype][jtype] != TWIST_NONE) {
magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall)
// omega_T (eq 29 of Marshall)
magtwist = relrot1*nx + relrot2*ny + relrot3*nz;
if (twist_model[itype][jtype] == TWIST_MARSHALL) {
k_twist = 0.5*k_tangential*a*a;; //eq 32
damp_twist = 0.5*damp_tangential*a*a;
@ -1525,11 +1623,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
damp_twist = twist_coeffs[itype][jtype][1];
mu_twist = twist_coeffs[itype][jtype][2];
}
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30)
// M_t torque (eq 30)
magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;
signtwist = (magtwist > 0) - (magtwist < 0);
Mtcrit = mu_twist*Fncrit;//critical torque (eq 44)
Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
if (fabs(magtortwist) > Mtcrit) {
magtortwist = -Mtcrit * signtwist; //eq 34
magtortwist = -Mtcrit * signtwist; // eq 34
}
}
@ -1578,8 +1677,8 @@ void PairGranular::unpack_forward_comm(int n, int first, double *buf)
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairGranular::memory_usage()
{
@ -1588,25 +1687,27 @@ double PairGranular::memory_usage()
}
/* ----------------------------------------------------------------------
mixing of Young's modulus (E)
mixing of Young's modulus (E)
------------------------------------------------------------------------- */
double PairGranular::mix_stiffnessE(double Eii, double Ejj, double poisii, double poisjj)
double PairGranular::mix_stiffnessE(double Eii, double Ejj,
double poisii, double poisjj)
{
return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj);
}
/* ----------------------------------------------------------------------
mixing of shear modulus (G)
mixing of shear modulus (G)
------------------------------------------------------------------------ */
double PairGranular::mix_stiffnessG(double Eii, double Ejj, double poisii, double poisjj)
double PairGranular::mix_stiffnessG(double Eii, double Ejj,
double poisii, double poisjj)
{
return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj));
}
/* ----------------------------------------------------------------------
mixing of everything else
mixing of everything else
------------------------------------------------------------------------- */
double PairGranular::mix_geom(double valii, double valjj)
@ -1616,10 +1717,11 @@ double PairGranular::mix_geom(double valii, double valjj)
/* ----------------------------------------------------------------------
Compute pull-off distance (beyond contact) for a given radius and atom type
compute pull-off distance (beyond contact) for a given radius and atom type
------------------------------------------------------------------------- */
double PairGranular::pulloff_distance(double radi, double radj, int itype, int jtype)
double PairGranular::pulloff_distance(double radi, double radj,
int itype, int jtype)
{
double E, coh, a, Reff;
Reff = radi*radj/(radi+radj);
@ -1631,11 +1733,12 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j
}
/* ----------------------------------------------------------------------
Transfer history during fix/neigh/history exchange
Only needed if any history entries i-j are not just negative of j-i entries
transfer history during fix/neigh/history exchange
only needed if any history entries i-j are not just negative of j-i entries
------------------------------------------------------------------------- */
void PairGranular::transfer_history(double* source, double* target) {
void PairGranular::transfer_history(double* source, double* target)
{
for (int i = 0; i < size_history; i++)
target[i] = history_transfer_factors[i]*source[i];
}

12
src/GRANULAR/pair_granular.h
View File

@ -66,29 +66,29 @@ class PairGranular : public Pair {
int size_history;
int *history_transfer_factors;
//Model choices
// model choices
int **normal_model, **damping_model;
int **tangential_model, **roll_model, **twist_model;
//History flags
// history flags
int normal_history, tangential_history, roll_history, twist_history;
//Indices of history entries
// indices of history entries
int normal_history_index;
int tangential_history_index;
int roll_history_index;
int twist_history_index;
//Per-type material coefficients
// per-type material coefficients
double **Emod, **poiss, **Gmod;
//Per-type coefficients, set in pair coeff command
// 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
// optional user-specified global cutoff, per-type user-specified cutoffs
double **cutoff_type;
double cutoff_global;

6
src/fix_neigh_history.cpp
View File

@ -408,7 +408,8 @@ void FixNeighHistory::pre_exchange_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues);
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues,jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}
@ -521,7 +522,8 @@ void FixNeighHistory::pre_exchange_no_newton()
m = npartner[j]++;
partner[j][m] = tag[i];
jvalues = &valuepartner[j][dnum*m];
if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues);
if (pair->nondefault_history_transfer)
pair->transfer_history(onevalues, jvalues);
else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
}
}

2
src/pair.h
View File

@ -98,7 +98,7 @@ class Pair : protected Pointers {
enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options
int beyond_contact, nondefault_history_transfer; //for granular styles
int beyond_contact, nondefault_history_transfer; // for granular styles
// KOKKOS host/device flag and data masks