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

Fixes several bugs in fix wall/gran, wall/gran/region:

1. Radius of curvature for curved regions was incorrectly used to compute wall-particle overlaps
2. Uninitialized values of rolling and twisting history could produce crashes in
cases that don't initialize these to 0 by default. These are now initialized to 0.
3. Fixed a bug with the use of 'NULL' for specification of the tangential stiffness for
wall/gran and wall/gran/region
parent 1e07ef6f
Loading
Loading
Loading
Loading
+0 −75
Original line number Diff line number Diff line
# pour cohesive particles on flat wall, let them flow
 
log		pour_wall_granular.log

atom_style	sphere
units		lj

###############################################
# Geometry-related parameters
###############################################

variable	boxx equal 50
variable	boxy equal 50
variable	boxz equal 20

variable	xc equal 0.5*${boxx}
variable	yc equal 0.5*${boxy}

###############################################
# Particle-related parameters
###############################################
variable	rlo equal 0.25
variable	rhi equal 0.5
variable	dlo equal 2.0*${rlo}
variable	dhi equal 2.0*${rhi}

variable	dens equal 1.0

variable skin equal 0.3*${rhi}

#############

region 		boxreg block 0 ${boxx} 0 ${boxy} 0 ${boxz}
create_box	2 boxreg
change_box	all boundary p p f

comm_modify 	vel yes

region		insreg cylinder z ${xc} ${yc} 5 15 ${boxz}

fix		1 all nve/sphere
fix		grav all gravity 10 vector 0 0 -1
fix		ins1 all pour 1000 1 3123 region insreg diam range ${dlo} ${dhi} dens ${dens} ${dens}
fix		ins2 all pour 1000 2 3123 region insreg diam range ${dlo} ${dhi} dens ${dens} ${dens}

comm_modify	vel yes

neighbor	${skin} bin
neigh_modify	delay 0 every 1 check yes

#pair_style	gran/hertz/history 1e5 1e5 1e3 1e3 0.5 1
#pair_coeff	* *

pair_style	granular 
pair_coeff	1 * hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5
#pair_coeff	2 2 jkr 1e5 1e4 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall
pair_coeff	2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
#pair_coeff	* * hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5

#fix		3 all wall/gran granular hertz/material 1e5 1e4 0.3 tangential mindlin 1e5 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall zplane 0 NULL 
fix		3 all wall/gran granular hertz/material 1e5 0.5 0.3 tangential mindlin 1e5 1.0 0.5 damping tsuji zplane 0 NULL 

thermo_style	custom step cpu atoms ke
thermo_modify	lost warn
thermo		100

timestep	0.001

dump		1 all custom 100 pour_wall_granular.dump id type radius mass x y z 

run		150000



	
+4 −6
Original line number Diff line number Diff line
@@ -122,6 +122,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
    iarg = 4;
    damping_model = VISCOELASTIC;
    roll_model = twist_model = NONE;
    tangential_history = roll_history = twist_history = 0;
    while (iarg < narg) {
      if (strcmp(arg[iarg], "hooke") == 0) {
        if (iarg + 2 >= narg)
@@ -226,7 +227,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
                         "stiffness requires a normal contact model "
                         "that specifies material properties");
            }
            tangential_coeffs[0] = 4*(2-poiss)*(1+poiss)/Emod;
            tangential_coeffs[0] = Emod/4*(2-poiss)*(1+poiss);
          } else {
            tangential_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt
          }
@@ -1061,7 +1062,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
                           double *contact)
{
  double fx,fy,fz,nx,ny,nz;
  double radsum,r,rinv;
  double r,rinv;
  double Reff, delta, dR, dR2;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -1095,11 +1096,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
  double shrmag,rsht;

  r = sqrt(rsq);
  radsum = rwall + radius;

  E = normal_coeffs[0];

  radsum = radius + rwall;
  if (rwall == 0) Reff = radius;
  else Reff = radius*rwall/(radius+rwall);

@@ -1122,7 +1120,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
  vn2 = ny*vnnr;
  vn3 = nz*vnnr;

  delta = radsum - r;
  delta = radius - r;
  dR = delta*Reff;
  if (normal_model == JKR) {
    history[0] = 1.0;
+2 −2

File changed.

Contains only whitespace changes.