Variable radius of region/sphere in LIGGGHTS

Submitted by Viraj on Tue, 09/15/2020 - 18:48


I have a really quick question. I have a spherical region, and I want its radius to change with time. On this
link , you can find the following statement:

"The radius value for style sphere and cylinder can be specified as an equal-style variable. If the value is a variable, it should be specified as v_name, where name is the variable name. In this case, the variable will be evaluated each timestep, and its value used to determine the radius of the region."

I am doing exactly the same thing (Input file can be found at the bottom). I could use either of the
following two commands to define the variable but I am getting errors (also mentioned):

variable b equal 0.1-(step/v_nsteps)*(0.1-0.09)
ERROR: Variable evaluation before simulation box is defined (../variable.cpp:1485)

variable b equal vdisplace(0.1,-0.01)
ERROR: Cannot use vdisplace in variable formula between runs (../variable.cpp:2783)

Can you please help me with defining a sphere region with a variable radius?

Thank you very much!

# declare basic details of the simulation
atom_style granular # the type of particles are granular
atom_modify map array # declare how atom id's look ups are performed
boundary p p p # boundaries in 3 directions are fixed (p = periodic)
newton off # not applying Newton's third law

communicate single vel yes # how interprocessor communication is done

units si # units are in SI

#variable b equal 0.1

#variable nsteps equal 1000000
variable rad equal 0.1-(step/v_nsteps)*(0.1-0.09)

variable b equal vdisplace(0.1,-0.01)

# declare what the simulation domain looks like by declaring a region
region domain sphere 0.0 0.0 0.0 v_rad units box
create_box 1 domain # create a box domain from the region

# declare details of cut-offs for neighbor lists, and how often they are rebuilt
neighbor 0.002 bin
neigh_modify delay 0

#Material properties required for new pair styles - can look up in documentation
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.9
fix m5 all property/global k_finnie peratomtypepair 1 1.0

#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * * # coefficients are given in prior command

timestep 0.00001 # set the time-step

fix gravi all gravity 0.0 vector 0.0 -1.0 0.0 # declare the #gravitational force

# set variable for going from inches to meters in stl files
# also how the stl files will be moved and scaled
variable inmeter equal 1
variable insmaller equal 0.001 #${inmeter}*1
variable impellerup equal 0.02032
variable inletup equal 0.2 #${inmeter}*0.2

# declare the pieces of the stl files
# the cylinder walls
# the spherical walls
fix sphere all mesh/surface file meshes/sphere.stl type 1 scale ${inmeter} #surface_vel -0.00001 spherical 0.0 0.0

# declare interaction between granules and the boundaries in the stl files
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes sphere

# distribution of particles for insertion
# smaller particles
fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2900 radius constant 0.005
# larger particles
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2900 radius constant 0.005
# ratio of large and small particles
fix pdd1 all particledistribution/discrete 32452843 2 pts1 0.5 pts2 0.5

# region and insertion
group nve_group region domain # declare groups of particles (all are in one group in this)

# particle insertion
# declare plane where particles will be inserted from
#fix inface all mesh/surface file meshes/InsertionFace.stl type 2 scale ${inmeter} move 0.0 ${inletup} 0.0
# declare how particles will be inserted
#fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
# nparticles 400 particlerate 10000 overlapcheck yes all_in no vel constant 0.0 -2.0 0.0 &
# insertion_face inface extrude_length 0.05

fix ins all insert/pack seed 123457 distributiontemplate pdd1 insert_every once overlapcheck yes volumefraction_region 0.6 region domain ntry_mc 10000

# apply nve (constant number, volume, energy) integration to all particles that are inserted as single particles (just a detail about ensembles)
fix integr nve_group nve/sphere

#fix 1 all deform 1 x vel -0.001 y vel -0.001 z vel -0.001

#output settings, include total thermal energy
compute 1 all erotate/sphere
compute str all stress/atom pair bond
compute myTemp all temp
compute pr all pressure myTemp
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no

#insert the first particles so that dump is not empty
run 1

#spin the impeller
#fix spin all move/mesh mesh impeller rotate origin 0. 0. 0. axis 0. 1. 0. period 1.

# output data
dump dmp all custom/vtk 2000 post/parts_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpcontainer all mesh/gran/VTK 2000 post/mesh_*.vtk vel area sphere
#dump dumpimpeller all mesh/gran/VTK 2000 post/mesh_impeller*.vtk vel area sphere

#insert particles
run 100000 upto
#run 1000 upto
unfix ins