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

Moved all model option syntax for pair granular to pair coeff command; 

added global cutoff option for pair style granular command; initial
write-up of documentation.
This commit is contained in:
dsbolin 2019-01-15 10:18:46 -07:00
parent 29dcdec875
commit ef803be08e
5 changed files with 538 additions and 595 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -0,0 +1,332 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
pair_style granular command :h3
pair_style granular/multi command :h3
[Syntax:]
pair_style style cutoff :pre
style = {granular} or {granular/multi} :ulb,l
cutoff = global cutoff (optional). See discussion below.
[Examples:]
pair_style granular
pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair coeff 2 2 hertz 200.0 20.0 tangential mindlin 300.0 50.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
pair_style granular/multi
pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair coeff 2 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall
pair coeff 1 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall
[Description:]
The {granular} styles support a variety of options for the normal, tangential, rolling and twisting
forces resulting from contact between two granular particles. The computed force depends on the combination
of choices for these models.
All model options 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
various combinations of particle types, as determined by the "pair_coeff"_pair_coeff.html command.
In the case of {granular}, coefficients
can vary between particle types, but model choices cannot. For instance, in the first
example above, the stiffness, damping, and tangential friction are different for type 1 - type 1 and type 2 - type 2 interactions, but
both 1-1 and 2-2 interactions must have the same model form, hence all keywords are identical between the two types. Cross-coefficients
for 1-2 interactions for the case of the {hertz} model above are set via simple geometric mixing rules. The {granular/multi}
style removes this restriction at a small cost in computational efficiency, so that different particle types
can potentially interact via different model forms. As shown in the second example,
1-1 interactions are based on a Hertzian contact model and 2-2 interactions are based on a {dmt} model (see below).
In the case that 1-1 and 2-2 interactions have different model forms, mixing of coefficients cannot be
determined, so 1-2 interactions must be explicitly defined via the pair coeff command, otherwise an error results.
The first required keyword for the pair coeff command is the normal contact model. Currently supported options and
the required arguments are:
{hooke} : k_n, damping
{hertz} : k_n, damping
{hertz/material} : E, damping, G
{dmt} : E, damping, G, cohesion
{jkr} : E, damping, G, cohesion
Here, k_n is spring stiffness, damping is a damping constant or a coefficient of restitution, depending on
the choice of damping model (see below), E and G are Young's modulus and shear modulus, in units of pressure,
and cohesion is a surface energy density, in units of energy/length^2.
For the {hooke} model, the normal component of force is given by:
:c,image(Eqs/hooke_normal.jpg)
For {hertz}, the normal force is given by:
:c,image{Eqs/hertz_normal.jpg}
For both [hooke] and [hertz], stiffness for unspecified cross-terms is given by simple geometric mixing
(e.g. if stiffness is specified for type 1 and type 2 particles as k_n_1 and k_n_2, respectively,
type 1 - type 2 contacts use a stiffness given by k_n_{12} = sqrt(k_n_1*k_n_2))
For {hertz/material}, the form is the same as above, but coefficients are computed differently, and mixing follows
a different rule based on shear modulus:
:c,image{Eqs/hertz_material_normal.jpg}
For {dmt}, the normal force is given by:
:c,image{Eqs/dmt_normal.jpg}
Where gamma is cohesion.
For {jkr}, the normal force is given by:
:c,image{Eqs/jkr_normal.jpg}
The same mixing rule for stiffness as for {hertz/material} is used by both the {dmt} and {jkr} models.
The tangential contact model must also be specified, which follows
the required {tangential} keyword. Currently supported options
and their required arguments are:
{no_history}: k_t, tangential_damping, friction coefficient
{mindlin}: k_t, tangential_damping, friction coefficient
For {no_history}, the tangential force is computed according to
The total force on a particle is the sum of the normal and tangential forces from all interactions. The tangential
force also induces a torque on both particles in a contacting pair. Additionally, rolling and twisting friction
models can also be applied, which may induce additional torques (but no force). The following options are
supported for the rolling friction model
The first required keyword
in the pair coeff command is the choice
of normal force contact model, for which current opitons are {hooke}, {hertz}
The {gran} styles use the following formulas for the frictional force
between two granular particles, as described in
"(Brilliantov)"_#Brilliantov, "(Silbert)"_#Silbert, and
"(Zhang)"_#Zhang3, when the distance r between two particles of radii
Ri and Rj is less than their contact distance d = Ri + Rj. There is
no force between the particles when r > d.
The two Hookean styles use this formula:
:c,image(Eqs/pair_gran_hooke.jpg)
The Hertzian style uses this formula:
:c,image(Eqs/pair_gran_hertz.jpg)
In both equations the first parenthesized term is the normal force
between the two particles and the second parenthesized term is the
tangential force. The normal force has 2 terms, a contact force and a
damping force. The tangential force also has 2 terms: a shear force
and a damping force. The shear force is a "history" effect that
accounts for the tangential displacement between the particles for the
duration of the time they are in contact. This term is included in
pair styles {hooke/history} and {hertz/history}, but is not included
in pair style {hooke}. The tangential damping force term is included
in all three pair styles if {dampflag} is set to 1; it is not included
if {dampflag} is set to 0.
The other quantities in the equations are as follows:
delta = d - r = overlap distance of 2 particles
Kn = elastic constant for normal contact
Kt = elastic constant for tangential contact
gamma_n = viscoelastic damping constant for normal contact
gamma_t = viscoelastic damping constant for tangential contact
m_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj
Delta St = tangential displacement vector between 2 particles \
which is truncated to satisfy a frictional yield criterion
n_ij = unit vector along the line connecting the centers of the 2 particles
Vn = normal component of the relative velocity of the 2 particles
Vt = tangential component of the relative velocity of the 2 particles :ul
The Kn, Kt, gamma_n, and gamma_t coefficients are specified as
parameters to the pair_style command. If a NULL is used for Kt, then
a default value is used where Kt = 2/7 Kn. If a NULL is used for
gamma_t, then a default value is used where gamma_t = 1/2 gamma_n.
The interpretation and units for these 4 coefficients are different in
the Hookean versus Hertzian equations.
The Hookean model is one where the normal push-back force for two
overlapping particles is a linear function of the overlap distance.
Thus the specified Kn is in units of (force/distance). Note that this
push-back force is independent of absolute particle size (in the
monodisperse case) and of the relative sizes of the two particles (in
the polydisperse case). This model also applies to the other terms in
the force equation so that the specified gamma_n is in units of
(1/time), Kt is in units of (force/distance), and gamma_t is in units
of (1/time).
The Hertzian model is one where the normal push-back force for two
overlapping particles is proportional to the area of overlap of the
two particles, and is thus a non-linear function of overlap distance.
Thus Kn has units of force per area and is thus specified in units of
(pressure). The effects of absolute particle size (monodispersity)
and relative size (polydispersity) are captured in the radii-dependent
pre-factors. When these pre-factors are carried through to the other
terms in the force equation it means that the specified gamma_n is in
units of (1/(time*distance)), Kt is in units of (pressure), and
gamma_t is in units of (1/(time*distance)).
Note that in the Hookean case, Kn can be thought of as a linear spring
constant with units of force/distance. In the Hertzian case, Kn is
like a non-linear spring constant with units of force/area or
pressure, and as shown in the "(Zhang)"_#Zhang3 paper, Kn = 4G /
(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E /
(2(1+nu)), and E = Young's modulus. Similarly, Kt = 4G / (2-nu).
(NOTE: in an earlier version of the manual, we incorrectly stated that
Kt = 8G / (2-nu).)
Thus in the Hertzian case Kn and Kt can be set to values that
corresponds to properties of the material being modeled. This is also
true in the Hookean case, except that a spring constant must be chosen
that is appropriate for the absolute size of particles in the model.
Since relative particle sizes are not accounted for, the Hookean
styles may not be a suitable model for polydisperse systems.
NOTE: In versions of LAMMPS before 9Jan09, the equation for Hertzian
interactions did not include the sqrt(RiRj/Ri+Rj) term and thus was
not as accurate for polydisperse systems. For monodisperse systems,
sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales all 4
coefficients: Kn, Kt, gamma_n, gamma_t. Thus you can set the values
of these 4 coefficients appropriately in the current code to reproduce
the results of a previous Hertzian monodisperse calculation. For
example, for the common case of a monodisperse system with particles
of diameter 1, all 4 of these coefficients should now be set 2x larger
than they were previously.
Xmu is also specified in the pair_style command and is the upper limit
of the tangential force through the Coulomb criterion Ft = xmu*Fn,
where Ft and Fn are the total tangential and normal force components
in the formulas above. Thus in the Hookean case, the tangential force
between 2 particles grows according to a tangential spring and
dash-pot model until Ft/Fn = xmu and is then held at Ft = Fn*xmu until
the particles lose contact. In the Hertzian case, a similar analogy
holds, though the spring is no longer linear.
NOTE: Normally, xmu should be specified as a fractional value between
0.0 and 1.0, however LAMMPS allows large values (up to 1.0e4) to allow
for modeling of systems which can sustain very large tangential
forces.
The effective mass {m_eff} is given by the formula above for two
isolated particles. If either particle is part of a rigid body, its
mass is replaced by the mass of the rigid body in the formula above.
This is determined by searching for a "fix rigid"_fix_rigid.html
command (or its variants).
For granular styles there are no additional coefficients to set for
each pair of atom types via the "pair_coeff"_pair_coeff.html command.
All settings are global and are made via the pair_style command.
However you must still use the "pair_coeff"_pair_coeff.html for all
pairs of granular atom types. For example the command
pair_coeff * * :pre
should be used if all atoms in the simulation interact via a granular
potential (i.e. one of the pair styles above is used). If a granular
potential is used as a sub-style of "pair_style
hybrid"_pair_hybrid.html, then specific atom types can be used in the
pair_coeff command to determine which atoms interact via a granular
potential.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
The "pair_modify"_pair_modify.html mix, shift, table, and tail options
are not relevant for granular pair styles.
These pair styles write their information to "binary restart
files"_restart.html, so a pair_style command does not need to be
specified in an input script that reads a restart file.
These pair styles can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. They do not support the
{inner}, {middle}, {outer} keywords.
The single() function of these pair styles returns 0.0 for the energy
of a pairwise interaction, since energy is not conserved in these
dissipative potentials. It also returns only the normal component of
the pairwise interaction force. However, the single() function also
calculates 10 extra pairwise quantities. The first 3 are the
components of the tangential force between particles I and J, acting
on particle I. The 4th is the magnitude of this tangential force.
The next 3 (5-7) are the components of the relative velocity in the
normal direction (along the line joining the 2 sphere centers). The
last 3 (8-10) the components of the relative velocity in the
tangential direction.
These extra quantities can be accessed by the "compute
pair/local"_compute_pair_local.html command, as {p1}, {p2}, ...,
{p10}.
:line
[Restrictions:]
All the granular pair styles are part of the GRANULAR package. It is
only enabled if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
These pair styles require that atoms store torque and angular velocity
(omega) as defined by the "atom_style"_atom_style.html. They also
require a per-particle radius is stored. The {sphere} atom style does
all of this.
This pair style requires you to use the "comm_modify vel
yes"_comm_modify.html command so that velocities are stored by ghost
atoms.
These pair styles will not restart exactly when using the
"read_restart"_read_restart.html command, though they should provide
statistically similar results. This is because the forces they
compute depend on atom velocities. See the
"read_restart"_read_restart.html command for more details.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Brilliantov)
[(Brilliantov)] Brilliantov, Spahn, Hertzsch, Poschel, Phys Rev E, 53,
p 5382-5392 (1996).
:link(Silbert)
[(Silbert)] Silbert, Ertas, Grest, Halsey, Levine, Plimpton, Phys Rev
E, 64, p 051302 (2001).
:link(Zhang3)
[(Zhang)] Zhang and Makse, Phys Rev E, 72, p 011301 (2005).

392
src/GRANULAR/pair_granular.cpp
View File

@ -1360,205 +1360,15 @@ void PairGranular::allocate()
void PairGranular::settings(int narg, char **arg)
{
if (narg < 5) error->all(FLERR,"Illegal pair_style command");
int iarg = 0;
//Some defaults
normal = HERTZ;
damping = VISCOELASTIC;
tangential = TANGENTIAL_MINDLIN;
roll = ROLL_NONE;
twist = TWIST_NONE;
tangential_history = 1;
if (narg == 1){
cutoff_global = force->numeric(FLERR,arg[0])
}
else{
cutoff_global = -1; //Will be set based on particle sizes, model choice
}
tangential_history = 0;
roll_history = twist_history = 0;
int normal_set, tangential_set;
normal_set = tangential_set = 0;
while (iarg < narg){
if (strcmp(arg[iarg], "hooke") == 0){
if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option");
normal = HOOKE;
memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_global[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_style command, not enough parameters provided for Hertz option");
normal = HERTZ;
normal_set = 1;
memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_global[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_style command, not enough parameters provided for Hertz option");
normal = HERTZ_MATERIAL;
normal_set = 1;
memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus)
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G (shear modulus)
iarg += num_coeffs+1;
}
else if (strcmp(arg[iarg], "dmt") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option");
normal = DMT;
normal_set = 1;
memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G
normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion
iarg += 5;
}
else if (strcmp(arg[iarg], "jkr") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option");
beyond_contact = 1;
normal = JKR;
normal_set = 1;
memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G
normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
iarg += 5;
}
else if (strcmp(arg[iarg], "damp_velocity") == 0){
damping = VELOCITY;
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){
damping = VISCOELASTIC;
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_tsuji") == 0){
damping = TSUJI;
iarg += 1;
}
else if (strstr(arg[iarg], "tangential") != NULL){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
tangential = TANGENTIAL_NOHISTORY;
tangential_set = 1;
}
else if (strstr(arg[iarg+1], "mindlin") != NULL){
tangential = TANGENTIAL_MINDLIN;
tangential_history = 1;
tangential_set = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized sliding friction model");
}
memory->create(tangential_coeffs_global, 3, "pair:tangential_coeffs_global");
tangential_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kt
tangential_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
tangential_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
else if (strstr(arg[iarg], "roll") != NULL){
if (strstr(arg[iarg+1], "none") != NULL){
roll = ROLL_NONE;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
roll = ROLL_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
roll = ROLL_SDS;
roll_history = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized rolling friction model");
}
memory->create(roll_coeffs_global, 3, "pair:roll_coeffs_global");
roll_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kR
roll_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR
roll_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
}
else if (strstr(arg[iarg], "twist") != NULL){
if (strstr(arg[iarg+1], "none") != NULL){
twist = TWIST_NONE;
iarg += 2;
}
else if (strstr(arg[iarg+1], "marshall") != NULL){
twist = TWIST_MARSHALL;
twist_history = 1;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model");
memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later
if (strstr(arg[iarg+1], "nohistory") != NULL){
twist = TWIST_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
twist = TWIST_SDS;
twist_history = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized twisting friction model");
}
memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global");
twist_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //ktwist
twist_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammatwist
twist_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
}
else error->all(FLERR, "Illegal pair_style granular command");
}
//Set all i-i entries, which may be replaced by pair coeff commands
//It may also make sense to consider removing all of the above, and only
// having the option for pair_coeff to set the parameters, similar to most LAMMPS pair styles
// The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax,
// where coeffs are set in the pair_style command, and a pair_coeff * * command is issued.
allocate();
double damp;
if (damping == TSUJI){
double cor = normal_coeffs_global[1];
damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+
27.467*pow(cor,4)-18.022*pow(cor,5)+
4.8218*pow(cor,6);
}
else damp = normal_coeffs_global[1];
for (int i = 1; i <= atom->ntypes; i++){
if (normal_set){
normal_coeffs[i][i][0] = normal_coeffs_global[0];
normal_coeffs[i][i][1] = damp;
if (normal != HOOKE && normal != HERTZ){
normal_coeffs[i][i][2] = normal_coeffs_global[2];
}
if ((normal == JKR) || (normal == DMT))
normal_coeffs[i][i][3] = normal_coeffs_global[3];
}
if(tangential_set){
for (int k = 0; k < 3; k++)
tangential_coeffs[i][i][k] = tangential_coeffs_global[k];
}
if (roll != ROLL_NONE)
for (int k = 0; k < 3; k++)
roll_coeffs[i][i][k] = roll_coeffs_global[k];
if (twist != TWIST_NONE && twist != TWIST_MARSHALL)
for (int k = 0; k < 3; k++)
twist_coeffs[i][i][k] = twist_coeffs_global[k];
if (normal_set && tangential_set) setflag[i][i] = 1;
}
normal_set = tangential_set = damping_set = roll_set = twist_set = 0;
}
/* ----------------------------------------------------------------------
@ -1567,12 +1377,6 @@ void PairGranular::settings(int narg, char **arg)
void PairGranular::coeff(int narg, char **arg)
{
int normal_set, damping_set, tangential_set, roll_set, twist_set;
double *normal_coeffs_local;
double *tangential_coeffs_local;
double *roll_coeffs_local;
double *twist_coeffs_local;
normal_coeffs_local = new double[4];
tangential_coeffs_local = new double[4];
roll_coeffs_local = new double[4];
@ -1587,13 +1391,12 @@ void PairGranular::coeff(int narg, char **arg)
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
normal_set = damping_set = tangential_set = roll_set = twist_set = 0;
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 (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent");
if (!normal_set) normal = HOOKE;
else if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types");
normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_set = 1;
@ -1602,7 +1405,8 @@ void PairGranular::coeff(int narg, char **arg)
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 (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent");
if (!normal_set) normal = HERTZ;
else if (normal_set && normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types");
normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_set = 1;
@ -1611,7 +1415,8 @@ void PairGranular::coeff(int narg, char **arg)
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 option");
if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent");
if (!normal_set) normal = HERTZ/MATERIAL;
else if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types");
normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E
normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G
@ -1620,7 +1425,8 @@ void PairGranular::coeff(int narg, char **arg)
}
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 (normal != DMT) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent");
if (!normal_set) normal = DMT;
else if (normal != DMT) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types");
normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E
normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G
@ -1631,7 +1437,8 @@ void PairGranular::coeff(int narg, char **arg)
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");
beyond_contact = 1;
if (normal != JKR) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent");
if (!normal_set) normal = JKR;
else if (normal != JKR) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types");
normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G
@ -1639,25 +1446,32 @@ void PairGranular::coeff(int narg, char **arg)
normal_set = 1;
iarg += 5;
}
else if (strcmp(arg[iarg], "damp_velocity") == 0){
if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent");
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){
if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent");
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_tsuji") == 0){
if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent");
iarg += 1;
}
else if (strstr(arg[iarg], "tangential") != NULL){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent");
else if (strcmp(arg[iarg], "damp") == 0){
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){
if (!damping_set) damping = VELOCITY;
else if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types");
}
else if (strstr(arg[iarg+1], "mindlin") != NULL){
if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent");;
else if (strcmp(arg[iarg+1], "viscoelastic") == 0){
if (!damping_set) damping = VISCOELASTIC;
else if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types");
}
else if (strcmp(arg[iarg+1], "tsuji") == 0){
if (!damping_set) damping = TSUJI;
if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types");
}
damping_set = 1;
iarg += 1;
}
else if (strcmp(arg[iarg], "tangential") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (strcmp(arg[iarg+1], "nohistory") == 0){
if (!tangential_set) tangential = TANGENTIAL_NOHISTORY;
else if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types");
}
else if (strcmp(arg[iarg+1], "mindlin") == 0){
if (!tangential_set) tangential = TANGENTIAL_MINDLIN;
else if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types");;
tangential_history = 1;
}
else{
@ -1669,19 +1483,22 @@ void PairGranular::coeff(int narg, char **arg)
tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
else if (strstr(arg[iarg], "rolling") != NULL){
else if (strcmp(arg[iarg], "rolling") == 0){
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strstr(arg[iarg+1], "none") != NULL){
if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent");
if (strcmp(arg[iarg+1], "none") == 0){
if (!roll_set) roll = ROLL_NONE;
else if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types");
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent");
if (strcmp(arg[iarg+1], "nohistory") == 0){
if (!roll_set) roll = ROLL_NOHISTORY;
else if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types");
}
else if (strstr(arg[iarg+1], "sds") != NULL){
if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent");
else if (strcmp(arg[iarg+1], "sds") == 0){
if (!roll_set) roll = ROLL_SDS;
else if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types");
roll_history = 1;
}
else{
@ -1694,24 +1511,28 @@ void PairGranular::coeff(int narg, char **arg)
iarg += 5;
}
}
else if (strstr(arg[iarg], "twist") != NULL){
else if (strcmp(arg[iarg], "twisting") == 0){
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strstr(arg[iarg+1], "none") != NULL){
if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent");
if (strcmp(arg[iarg+1], "none") == 0){
if (!twist_set) twist = TWIST_NOHISTORY;
else if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types");
iarg += 2;
}
else if (strstr(arg[iarg+1], "marshall") != NULL){
if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent");
else if (strcmp(arg[iarg+1], "marshall") == 0){
if (!twist_set) twist = TWIST_MARSHALL;
else if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types");
twist_history = 1;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent");
if (!twist_set) twist = TWIST_NOHISTORY;
else if (strcmp(arg[iarg+1], "nohistory") == 0){
if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types");
}
else if (strstr(arg[iarg+1], "sds") != NULL){
if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent");
else if (strcmp(arg[iarg+1], "sds") == 0){
if (!twist_set) twist = TWIST_SDS;
else if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types");
twist_history = 1;
}
else{
@ -1727,6 +1548,15 @@ void PairGranular::coeff(int narg, char **arg)
else error->all(FLERR, "Illegal pair coeff command");
}
//It is an error not to specify normal or tangential model
if (!normal_set) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model");
if (!tangential_set) error->all(FLERR, "Illegal pair coeff command, must specify tangential contact model");
//If unspecified, set damping to VISCOELASTIC, twist/roll to NONE (cannot be changed by subsequent pair_coeff commands)
if (!damping_set) damping = VISCOELASTIC;
if (!roll_set) roll = ROLL_NONE;
if (!twist_set) twist = TWIST_NONE;
int count = 0;
double damp;
if (damping == TSUJI){
@ -1740,37 +1570,28 @@ void PairGranular::coeff(int narg, char **arg)
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
if (normal_set){
normal_coeffs[i][j][0] = normal_coeffs_local[0];
normal_coeffs[i][j][1] = damp;
if (normal != HERTZ && normal != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2];
if ((normal == JKR) || (normal == DMT))
normal_coeffs[i][j][3] = normal_coeffs_local[3];
}
if (tangential_set){
normal_coeffs[i][j][0] = normal_coeffs_local[0];
normal_coeffs[i][j][1] = damp;
if (normal != HERTZ && normal != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2];
if ((normal == JKR) || (normal == DMT))
normal_coeffs[i][j][3] = normal_coeffs_local[3];
for (int k = 0; k < 3; k++)
tangential_coeffs[i][j][k] = tangential_coeffs_local[k];
if (roll != ROLL_NONE)
for (int k = 0; k < 3; k++)
tangential_coeffs[i][j][k] = tangential_coeffs_local[k];
}
if (roll_set){
if (roll != ROLL_NONE)
for (int k = 0; k < 3; k++)
roll_coeffs[i][j][k] = roll_coeffs_local[k];
}
if (twist_set){
if (twist != TWIST_NONE && twist != TWIST_MARSHALL)
for (int k = 0; k < 3; k++)
twist_coeffs[i][j][k] = twist_coeffs_local[k];
}
roll_coeffs[i][j][k] = roll_coeffs_local[k];
if (twist != TWIST_NONE && twist != TWIST_MARSHALL)
for (int k = 0; k < 3; k++)
twist_coeffs[i][j][k] = twist_coeffs_local[k];
setflag[i][j] = 1;
count++;
}
}
delete[] normal_coeffs_local;
delete[] tangential_coeffs_local;
delete[] roll_coeffs_local;
delete[] twist_coeffs_local;
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
@ -1959,20 +1780,25 @@ double PairGranular::init_one(int i, int j)
// 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.
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
cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
}
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]);
cutmax = MAX(cutmax,2.0*maxrad_frozen[k]);
if (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
cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
}
cutoff = cutmax;
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]);
cutmax = MAX(cutmax,2.0*maxrad_frozen[k]);
}
cutoff = cutmax;
}
}
else{
cutoff = cutoff_global;
}
return cutoff;
}

13
src/GRANULAR/pair_granular.h
View File

@ -86,19 +86,18 @@ private:
//Indices of history entries
int tangential_history_index, roll_history_index, twist_history_index;
//Coefficients declared in pair style command, used as default unless
// overwritten in pair coeff command
double *normal_coeffs_global;
double *tangential_coeffs_global;
double *roll_coeffs_global;
double *twist_coeffs_global;
//Flags for whether model choices have been set
int normal_set, tangential_set, damping_set, roll_set, twist_set;
//Per-type coefficients declared 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
double global_cutoff;
double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj);
double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj);
double mix_geom(double valii, double valjj);

379
src/GRANULAR/pair_granular_multi.cpp
View File

@ -654,217 +654,15 @@ void PairGranularMulti::allocate()
void PairGranularMulti::settings(int narg, char **arg)
{
if (narg < 5) error->all(FLERR,"Illegal pair_style command");
if (narg == 1){
cutoff_global = force->numeric(FLERR,arg[0])
}
else{
cutoff_global = -1; //Will be set based on particle sizes, model choice
}
int iarg = 0;
//Some defaults
normal_global = HERTZ;
damping_global = VISCOELASTIC;
tangential_global = TANGENTIAL_MINDLIN;
roll_global = ROLL_NONE;
twist_global = TWIST_NONE;
tangential_history = 1;
tangential_history = 0;
roll_history = twist_history = 0;
int normal_set, tangential_set;
normal_set = tangential_set = 0;
while (iarg < narg){
if (strcmp(arg[iarg], "hooke") == 0){
if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option");
normal_global = HOOKE;
normal_set = 1;
memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_global[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_style command, not enough parameters provided for Hertz option");
normal_global = HERTZ;
normal_set = 1;
memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn
normal_coeffs_global[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_style command, not enough parameters provided for Hertz option");
normal_global = HERTZ_MATERIAL;
normal_set = 1;
memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus)
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G (shear modulus)
iarg += num_coeffs+1;
}
else if (strcmp(arg[iarg], "dmt") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option");
normal_global = DMT;
normal_set = 1;
memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G
normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion
iarg += 5;
}
else if (strcmp(arg[iarg], "jkr") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option");
beyond_contact = 1;
normal_global = JKR;
normal_set = 1;
memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global");
normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E
normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping
normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G
normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
iarg += 5;
}
else if (strcmp(arg[iarg], "damp_velocity") == 0){
damping_global = VELOCITY;
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){
damping_global = VISCOELASTIC;
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_tsuji") == 0){
damping_global = TSUJI;
iarg += 1;
}
else if (strstr(arg[iarg], "tangential") != NULL){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
tangential_global = TANGENTIAL_NOHISTORY;
tangential_set = 1;
}
else if (strstr(arg[iarg+1], "mindlin") != NULL){
tangential_global = TANGENTIAL_MINDLIN;
tangential_history = 1;
tangential_set = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized sliding friction model");
}
memory->create(tangential_coeffs_global, 3, "pair:tangential_coeffs_global");
tangential_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kt
tangential_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
tangential_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
else if (strstr(arg[iarg], "roll") != NULL){
if (strstr(arg[iarg+1], "none") != NULL){
roll_global = ROLL_NONE;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
roll_global = ROLL_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
roll_global = ROLL_SDS;
roll_history = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized rolling friction model");
}
memory->create(roll_coeffs_global, 3, "pair:roll_coeffs_global");
roll_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kR
roll_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR
roll_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
}
else if (strstr(arg[iarg], "twist") != NULL){
if (strstr(arg[iarg+1], "none") != NULL){
twist_global = TWIST_NONE;
iarg += 2;
}
else if (strstr(arg[iarg+1], "marshall") != NULL){
twist_global = TWIST_MARSHALL;
twist_history = 1;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model");
memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later
if (strstr(arg[iarg+1], "nohistory") != NULL){
twist_global = TWIST_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
twist_global = TWIST_SDS;
twist_history = 1;
}
else{
error->all(FLERR, "Illegal pair_style command, unrecognized twisting friction model");
}
memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global");
twist_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //ktwist
twist_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammatwist
twist_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
}
else error->all(FLERR, "Illegal pair_style granular command");
}
//Set all i-i entries, which may be replaced by pair coeff commands
//It may also make sense to consider removing all of the above, and only
// having the option for pair_coeff to set the parameters, similar to most LAMMPS pair styles
// The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax,
// where coeffs are set in the pair_style command, and a pair_coeff * * command is issued.
//Other option is to have two pair styles, e.g. pair granular and pair granular/multi,
// where granular/multi allows per-type coefficients, pair granular does not (this would also
// allow minor speed-up by templating pair granular)
allocate();
double damp;
for (int i = 1; i <= atom->ntypes; i++){
normal[i][i] = normal_global;
damping[i][i] = damping_global;
tangential[i][i] = tangential_global;
roll[i][i] = roll_global;
twist[i][i] = twist_global;
if (normal_set){
if (damping_global == TSUJI){
double cor = normal_coeffs_global[1];
damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+
27.467*pow(cor,4)-18.022*pow(cor,5)+
4.8218*pow(cor,6);
}
else damp = normal_coeffs_global[1];
normal_coeffs[i][i][0] = normal_coeffs_global[0];
normal_coeffs[i][i][1] = damp;
if (normal[i][i] != HOOKE && normal[i][i] != HERTZ){
normal_coeffs[i][i][2] = normal_coeffs_global[2];
}
if ((normal_global == JKR) || (normal_global == DMT))
normal_coeffs[i][i][3] = normal_coeffs_global[3];
}
if(tangential_set){
tangential[i][i] = tangential_global;
for (int k = 0; k < 3; k++)
tangential_coeffs[i][i][k] = tangential_coeffs_global[k];
}
roll[i][i] = roll_global;
if (roll_global != ROLL_NONE)
for (int k = 0; k < 3; k++)
roll_coeffs[i][i][k] = roll_coeffs_global[k];
twist[i][i] = twist_global;
if (twist_global != TWIST_NONE && twist_global != TWIST_MARSHALL)
for (int k = 0; k < 3; k++)
twist_coeffs[i][i][k] = twist_coeffs_global[k];
if (normal_set && tangential_set) setflag[i][i] = 1;
}
}
/* ----------------------------------------------------------------------
@ -874,10 +672,6 @@ void PairGranularMulti::settings(int narg, char **arg)
void PairGranularMulti::coeff(int narg, char **arg)
{
int normal_local, damping_local, tangential_local, roll_local, twist_local;
double *normal_coeffs_local;
double *tangential_coeffs_local;
double *roll_coeffs_local;
double *twist_coeffs_local;
normal_coeffs_local = new double[4];
tangential_coeffs_local = new double[4];
@ -893,8 +687,10 @@ void PairGranularMulti::coeff(int narg, char **arg)
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
normal_local = tangential_local = roll_local = twist_local = -1;
damping_local = -1;
//Defaults
normal_local = tangential_local = -1;
roll_local = twist_local = 0;
damping_local = VISCOELASTIC;
int iarg = 2;
while (iarg < narg){
@ -941,24 +737,27 @@ void PairGranularMulti::coeff(int narg, char **arg)
normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
iarg += 5;
}
else if (strcmp(arg[iarg], "damp_velocity") == 0){
damping_local = VELOCITY;
iarg += 1;
else if (strcmp(arg[iarg], "damp") == 0){
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_local = VELOCITY;
iarg += 1;
}
else if (strcmp(arg[iarg+1], "viscoelastic") == 0){
damping_local = VISCOELASTIC;
iarg += 1;
}
else if (strcmp(arg[iarg], "tsuji") == 0){
damping_local = TSUJI;
iarg += 1;
}
}
else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){
damping_local = VISCOELASTIC;
iarg += 1;
}
else if (strcmp(arg[iarg], "damp_tsuji") == 0){
damping_local = TSUJI;
iarg += 1;
}
else if (strstr(arg[iarg], "tangential") != NULL){
else if (strcmp(arg[iarg], "tangential") == 0){
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (strcmp(arg[iarg+1], "nohistory") == 0){
tangential_local = TANGENTIAL_NOHISTORY;
}
else if (strstr(arg[iarg+1], "mindlin") != NULL){
else if (strcmp(arg[iarg+1], "mindlin") == 0){
tangential_local = TANGENTIAL_MINDLIN;
tangential_history = 1;
}
@ -970,18 +769,18 @@ void PairGranularMulti::coeff(int narg, char **arg)
tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
iarg += 5;
}
else if (strstr(arg[iarg], "rolling") != NULL){
else if (strcmp(arg[iarg], "rolling") == 0){
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strstr(arg[iarg+1], "none") != NULL){
if (strcmp(arg[iarg+1], "none") == 0){
roll_local = ROLL_NONE;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (strcmp(arg[iarg+1], "nohistory") == 0){
roll_local = ROLL_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
else if (strcmp(arg[iarg+1], "sds") == 0){
roll_local = ROLL_SDS;
roll_history = 1;
}
@ -994,23 +793,23 @@ void PairGranularMulti::coeff(int narg, char **arg)
iarg += 5;
}
}
else if (strstr(arg[iarg], "twist") != NULL){
else if (strcmp(arg[iarg], "twisting") == 0){
if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
if (strstr(arg[iarg+1], "none") != NULL){
if (strcmp(arg[iarg+1], "none") == 0){
twist_local = TWIST_NONE;
iarg += 2;
}
else if (strstr(arg[iarg+1], "marshall") != NULL){
else if (strcmp(arg[iarg+1], "marshall") == 0){
twist_local = TWIST_MARSHALL;
twist_history = 1;
iarg += 2;
}
else{
if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model");
if (strstr(arg[iarg+1], "nohistory") != NULL){
if (strcmp(arg[iarg+1], "nohistory") == 0){
twist_local = TWIST_NOHISTORY;
}
else if (strstr(arg[iarg+1], "sds") != NULL){
else if (strcmp(arg[iarg+1], "sds") == 0){
twist_local = TWIST_SDS;
twist_history = 1;
}
@ -1026,66 +825,49 @@ void PairGranularMulti::coeff(int narg, char **arg)
else error->all(FLERR, "Illegal pair coeff command");
}
//It is an error not to specify normal or tangential model
if ((normal_set < 0) || (tangential_set < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model");))
int count = 0;
double damp;
if (damping_local >= 0){
if (normal_local == -1)
error->all(FLERR, "Illegal pair_coeff command, must specify normal model when setting damping model");
if (damping_local == TSUJI){
double cor;
cor = normal_coeffs_local[1];
damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+
27.467*pow(cor,4)-18.022*pow(cor,5)+
4.8218*pow(cor,6);
}
else damp = normal_coeffs_local[1];
if (damping_local == TSUJI){
double cor;
cor = normal_coeffs_local[1];
damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+
27.467*pow(cor,4)-18.022*pow(cor,5)+
4.8218*pow(cor,6);
}
else damp = normal_coeffs_local[1];
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
if (normal_local >= 0){
normal[i][j] = normal_local;
normal_coeffs[i][j][0] = normal_coeffs_local[0];
if (damping_local == -1){
damp = normal_coeffs_global[1];
}
normal_coeffs[i][j][1] = damp;
if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2];
if ((normal_local == JKR) || (normal_local == DMT))
normal_coeffs[i][j][3] = normal_coeffs_local[3];
}
if (damping_local >= 0){
damping[i][j] = damping_local;
}
if (tangential_local >= 0){
tangential[i][j] = tangential_local;
normal[i][j] = normal_local;
normal_coeffs[i][j][0] = normal_coeffs_local[0];
normal_coeffs[i][j][1] = damp;
if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2];
if ((normal_local == JKR) || (normal_local == DMT))
normal_coeffs[i][j][3] = normal_coeffs_local[3];
damping[i][j] = damping_local;
tangential[i][j] = tangential_local;
for (int k = 0; k < 3; k++)
tangential_coeffs[i][j][k] = tangential_coeffs_local[k];
}
if (roll_local >= 0){
roll[i][j] = roll_local;
if (roll_local != ROLL_NONE)
for (int k = 0; k < 3; k++)
roll_coeffs[i][j][k] = roll_coeffs_local[k];
}
if (twist_local >= 0){
twist[i][j] = twist_local;
roll[i][j] = roll_local;
if (roll_local != ROLL_NONE)
for (int k = 0; k < 3; k++)
roll_coeffs[i][j][k] = roll_coeffs_local[k];
twist[i][j] = twist_local;
if (twist_local != TWIST_NONE && twist_local != TWIST_MARSHALL)
for (int k = 0; k < 3; k++)
twist_coeffs[i][j][k] = twist_coeffs_local[k];
}
if (normal_local >= 0 && tangential_local >= 0) setflag[i][j] = 1;
setflag[i][j] = 1;
count++;
}
}
delete[] normal_coeffs_local;
delete[] tangential_coeffs_local;
delete[] roll_coeffs_local;
delete[] twist_coeffs_local;
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
@ -1286,20 +1068,25 @@ double PairGranularMulti::init_one(int i, int j)
// 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.
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
cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
}
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]);
cutmax = MAX(cutmax,2.0*maxrad_frozen[k]);
if (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
cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
}
cutoff = cutmax;
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]);
cutmax = MAX(cutmax,2.0*maxrad_frozen[k]);
}
cutoff = cutmax;
}
}
else{
cutoff = cutoff_global;
}
return cutoff;
}

17
src/GRANULAR/pair_granular_multi.h
View File

@ -65,27 +65,26 @@ public:
private:
int size_history;
//Per-type models
//Models
int **normal, **damping, **tangential, **roll, **twist;
int normal_global, damping_global;
int tangential_global, roll_global, twist_global;
//History flags
int tangential_history, roll_history, twist_history;
//Indices of history entries
int tangential_history_index;
int roll_history_index;
int twist_history_index;
double *normal_coeffs_global;
double *tangential_coeffs_global;
double *roll_coeffs_global;
double *twist_coeffs_global;
//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
double cutoff_global;
double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj);
double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj);
double mix_geom(double valii, double valjj);