Commit 4e26ca29 authored by Dan S. Bolintineanu's avatar Dan S. Bolintineanu
Browse files

Changes to new generalized granular pair styles and fix wall/gran 

-Clean-up of unused variables in code
-Bug fix for single method of pair granular
-Changes to fix wall/gran to fix issues with JKR
-Doc page updates for fix wall/gran and fix wall/gran/region
parent ff795e76

doc/src/fix_wall_gran.txt

+34 −19
Original line number Diff line number Diff line
@@ -11,18 +11,21 @@ fix wall/gran/omp command :h3 


[Syntax:]



fix ID group-ID wall/gran fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle args keyword values ... :pre

fix ID group-ID wall/gran fstyle fstyle_params wallstyle args keyword values ... :pre



ID, group-ID are documented in "fix"_fix.html command :ulb,l

wall/gran = style name of this fix command :l

fstyle = style of force interactions between particles and wall :l

  possible choices: hooke, hooke/history, hertz/history :pre

Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l

Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l

gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l

gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l

xmu = static yield criterion (unitless value between 0.0 and 1.0e4) :l

dampflag = 0 or 1 if tangential damping force is excluded or included :l

  possible choices: hooke, hooke/history, hertz/history, granular :pre

fstyle_params = parameters associated with force interaction style :l

  For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are: 

	Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) 

	Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) 

	gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) 

	gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)  

	xmu = static yield criterion (unitless value between 0.0 and 1.0e4) 

	dampflag = 0 or 1 if tangential damping force is excluded or included :pre

  For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre

wallstyle = {xplane} or {yplane} or {zplane} or {zcylinder} :l

args = list of arguments for a particular style :l

  {xplane} or {yplane} or {zplane} args = lo hi
@@ -44,7 +47,10 @@ keyword = {wiggle} or {shear} :l 


fix 1 all wall/gran hooke  200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0

fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL

fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 :pre

fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 

fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 zplane 0.0 NULL

fix 4 all wall/gran granular jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall zcylinder 15.0 wiggle z 3.0 2.0

fix 5 all wall/gran granular dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall zplane 0.0 NULL :pre



[Description:]
@@ -54,31 +60,39 @@ close enough to touch it. 


The nature of the wall/particle interactions are determined by the

{fstyle} setting.  It can be any of the styles defined by the

"pair_style granular"_pair_gran.html commands.  Currently this is

{hooke}, {hooke/history}, or {hertz/history}.  The equation for the

"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 granular"_pair_gran.html doc

page, in the limit of one of the two particles going to infinite

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 just the radius of the

particle.

effective radius of contact = RiRj/Ri+Rj is set to the radius

of the particle.



The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}

have the same meaning and units as those specified with the

"pair_style granular"_pair_gran.html commands.  This means a NULL can

"pair_style gran/*"_pair_gran.html commands.  This means a NULL can

be used for either {Kt} or {gamma_t} as described on that page.  If a

NULL is used for {Kt}, then a default value is used where {Kt} = 2/7

{Kn}.  If a NULL is used for {gamma_t}, then a default value is used

where {gamma_t} = 1/2 {gamma_n}.



All the model choices for cohesion, tangential friction, rolling friction

and twisting friction supported by the "pair_style granular"_pair_granular.html

through its {pair_coeff} command are also supported for walls. These are discussed 

in greater detail on the doc page for "pair_style granular"_pair_granular.html.



Note that you can choose a different force styles and/or different

values for the 6 wall/particle coefficients than for particle/particle

values for the wall/particle coefficients than for particle/particle

interactions.  E.g. if you wish to model the wall as a different

material.



NOTE: As discussed on the doc page for "pair_style

granular"_pair_gran.html, versions of LAMMPS before 9Jan09 used a

gran/*"_pair_gran.html, versions of LAMMPS before 9Jan09 used a

different equation for Hertzian interactions.  This means Hertizian

wall/particle interactions have also changed.  They now include a

sqrt(radius) term which was not present before.  Also the previous
@@ -188,6 +202,7 @@ Any dimension (xyz) that has a granular wall must be non-periodic. 


"fix move"_fix_move.html,

"fix wall/gran/region"_fix_wall_gran_region.html,

"pair_style granular"_pair_gran.html

"pair_style gran/*"_pair_gran.html

"pair_style granular"_pair_granular.html



[Default:] none

doc/src/fix_wall_gran_region.txt

+29 −14
Original line number Diff line number Diff line
@@ -10,24 +10,30 @@ fix wall/gran/region command :h3 


[Syntax:]



fix ID group-ID wall/gran/region fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle regionID :pre

fix ID group-ID wall/gran/region fstyle fstyle_params wallstyle regionID :pre



ID, group-ID are documented in "fix"_fix.html command :ulb,l

wall/region = style name of this fix command :l

fstyle = style of force interactions between particles and wall :l

  possible choices: hooke, hooke/history, hertz/history :pre

Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l

Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l

gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l

gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l

xmu = static yield criterion (unitless value between 0.0 and 1.0e4) :l

dampflag = 0 or 1 if tangential damping force is excluded or included :l

  possible choices: hooke, hooke/history, hertz/history, granular :pre

fstyle_params = parameters associated with force interaction style :l

  For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are: 

	Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) 

	Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) 

	gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) 

	gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)  

	xmu = static yield criterion (unitless value between 0.0 and 1.0e4) 

	dampflag = 0 or 1 if tangential damping force is excluded or included :pre

  For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre

wallstyle = region (see "fix wall/gran"_fix_wall_gran.html for options for other kinds of walls) :l

region-ID = region whose boundary will act as wall :l,ule



[Examples:]



fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone :pre

fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone 

fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 region myBox

fix 4 all wall/gran/region granular jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone

fix 5 all wall/gran/region granular dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall region myCone :pre



[Description:]
@@ -122,11 +128,14 @@ to make the two faces differ by epsilon in their position. 


The nature of the wall/particle interactions are determined by the

{fstyle} setting.  It can be any of the styles defined by the

"pair_style granular"_pair_gran.html commands.  Currently this is

{hooke}, {hooke/history}, or {hertz/history}.  The equation for the

"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 granular"_pair_gran.html doc

page, but the effective radius is calculated using the radius of 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,
@@ -140,12 +149,18 @@ particle. 


The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}

have the same meaning and units as those specified with the

"pair_style granular"_pair_gran.html commands.  This means a NULL can

"pair_style gran/*"_pair_gran.html commands.  This means a NULL can

be used for either {Kt} or {gamma_t} as described on that page.  If a

NULL is used for {Kt}, then a default value is used where {Kt} = 2/7

{Kn}.  If a NULL is used for {gamma_t}, then a default value is used

where {gamma_t} = 1/2 {gamma_n}.





All the model choices for cohesion, tangential friction, rolling friction

and twisting friction supported by the "pair_style granular"_pair_granular.html

through its {pair_coeff} command are also supported for walls. These are discussed 

in greater detail on the doc page for "pair_style granular"_pair_granular.html.



Note that you can choose a different force styles and/or different

values for the 6 wall/particle coefficients than for particle/particle

interactions.  E.g. if you wish to model the wall as a different

src/GRANULAR/fix_wall_gran.cpp

+6 −10
Original line number Diff line number Diff line
@@ -1021,9 +1021,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, 
    double radius, double meff, double *history,

    double *contact)

{

  int i,j,ii,jj,inum,jnum,itype,jtype;

 double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz;

 double radi,radj,radsum,r,rinv,rsqinv;

 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;
@@ -1035,17 +1034,17 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, 
 double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;

 double fs, fs1, fs2, fs3;



 double mi,mj,damp,ccel,tor1,tor2,tor3;

 double tor1,tor2,tor3;

 double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv;



 //For JKR

 double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E;

 double R2, coh, F_pulloff, a, a2, E;

 double t0, t1, t2, t3, t4, t5, t6;

 double sqrt1, sqrt2, sqrt3, sqrt4;

 double sqrt1, sqrt2, sqrt3;



 //Rolling

 double k_roll, damp_roll;

 double roll1, roll2, roll3, torroll1, torroll2, torroll3;

 double torroll1, torroll2, torroll3;

 double rollmag, rolldotn, scalefac;

 double fr, fr1, fr2, fr3;
@@ -1055,9 +1054,6 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, 
 double tortwist1, tortwist2, tortwist3;



 double shrmag,rsht;

 int *ilist,*jlist,*numneigh,**firstneigh;

 int *touch,**firsttouch;

 double *allhistory,**firsthistory;



 r = sqrt(rsq);

 radsum = rwall + radius;

src/GRANULAR/pair_gran_hooke_history.cpp

+0 −4
Original line number Diff line number Diff line
@@ -349,10 +349,6 @@ void PairGranHookeHistory::settings(int narg, char **arg) 
{

  if (narg != 6) error->all(FLERR,"Illegal pair_style command");











  kn = force->numeric(FLERR,arg[0]);

  if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0;

  else kt = force->numeric(FLERR,arg[1]);

src/GRANULAR/pair_granular.cpp

+16 −34
Original line number Diff line number Diff line
@@ -130,7 +130,7 @@ void PairGranular::compute(int eflag, int vflag) 
{

  int i,j,ii,jj,inum,jnum,itype,jtype;

  double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz;

  double radi,radj,radsum,rsq,r,rinv,rsqinv;

  double radi,radj,radsum,rsq,r,rinv;

  double Reff, delta, dR, dR2;



  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -140,19 +140,19 @@ void PairGranular::compute(int eflag, int vflag) 
  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 fs, fs1, fs2, fs3, tor1, tor2, tor3;



  double mi,mj,meff,damp,ccel,tor1,tor2,tor3;

  double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv;

  double mi,mj,meff;

  double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;



  //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, sqrt4;

  double sqrt1, sqrt2, sqrt3;



  //Rolling

  double k_roll, damp_roll;

  double roll1, roll2, roll3, torroll1, torroll2, torroll3;

  double torroll1, torroll2, torroll3;

  double rollmag, rolldotn, scalefac;

  double fr, fr1, fr2, fr3;
@@ -204,7 +204,6 @@ void PairGranular::compute(int eflag, int vflag) 
  double *rmass = atom->rmass;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;
@@ -465,9 +464,6 @@ void PairGranular::compute(int eflag, int vflag) 
          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;

          vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3);

          if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag;

          else vrlmaginv = 0.0;



          int rhist0 = roll_history_index;

          int rhist1 = rhist0 + 1;
@@ -1192,12 +1188,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    double rsq, double factor_coul, double factor_lj, double &fforce)

{

  double radi,radj,radsum;

  double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, Reff;

  double r,rinv,delx,dely,delz, nx, ny, nz, Reff;

  double dR, dR2;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;

  double vtr1,vtr2,vtr3,vrel;

  double mi,mj,meff,damp,ccel,tor1,tor2,tor3;

  double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv;

  double mi,mj,meff;

  double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;



  double knfac, damp_normal, damp_normal_prefactor;

  double k_tangential, damp_tangential;
@@ -1207,25 +1203,22 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  //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, sqrt4;

  double sqrt1, sqrt2, sqrt3;





  //Rolling

  double k_roll, damp_roll;

  double roll1, roll2, roll3, torroll1, torroll2, torroll3;

  double rollmag, rolldotn, scalefac;

  double rollmag;

  double fr, fr1, fr2, fr3;



  //Twisting

  double k_twist, damp_twist, mu_twist;

  double signtwist, magtwist, magtortwist, Mtcrit;

  double tortwist1, tortwist2, tortwist3;



  double shrmag,rsht;

  double shrmag;

  int jnum;

  int *ilist,*jlist,*numneigh,**firstneigh;

  int *touch,**firsttouch;

  double *history,*allhistory,**firsthistory;

  int *jlist;

  double *history,*allhistory;



  double *radius = atom->radius;

  radi = radius[i];
@@ -1312,8 +1305,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  // if I or J part of rigid body, use body mass

  // if I or J is frozen, meff is other particle



  int *type = atom->type;



  mi = rmass[i];

  mj = rmass[j];

  if (fix_rigid) {
@@ -1348,11 +1339,11 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
  else{

    knfac = E;

    a = sqrt(dR);

    Fne = knfac*delta;

    if (normal_model[itype][jtype] != HOOKE){

      Fne *= a;

      knfac *= a;

    }

    Fne = knfac*delta;

    if (normal_model[itype][jtype] == DMT)

      Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;

  }
@@ -1471,9 +1462,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    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;

    vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3);

    if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag;

    else vrlmaginv = 0.0;



    int rhist0 = roll_history_index;

    int rhist1 = rhist0 + 1;
@@ -1484,8 +1472,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
        history[rhist1]*history[rhist1] +

        history[rhist2]*history[rhist2]);



    rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;



    k_roll = roll_coeffs[itype][jtype][0];

    damp_roll = roll_coeffs[itype][jtype][1];

    fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
@@ -1498,9 +1484,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    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;
@@ -1528,7 +1511,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, 
    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

    }

  }
@@ -1621,7 +1603,7 @@ double PairGranular::mix_geom(double valii, double valjj) 


double PairGranular::pulloff_distance(double radi, double radj, int itype, int jtype)

{

  double E, coh, a, delta_pulloff, Reff;

  double E, coh, a, Reff;

  Reff = radi*radj/(radi+radj);

  if (Reff <= 0) return 0;

  coh = normal_coeffs[itype][itype][3];

CRA Git | Maintained and supported by SUSTech CRA and CCSE