Hi all,
I am trying to do simulation on a rotary drum, simulating the rotation by changing the gravity force. I hope you could help me. Here is my code:
atom_style granular
atom_modify map array
boundary p p p
echo screen
newton off
units si
communicate single vel yes
region reg block 0 2.4 0 0.314 0 0.314 units box
create_box 1 reg
neighbor 0.02 bin
neigh_modify delay 0
fix m1 all property/global youngsModulus peratomtype 2.1e11
fix m2 all property/global poissonsRatio peratomtype 0.27
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.95
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
pair_style gran/hertz/history
pair_coeff * *
timestep 0.00001
variable g equal 0.00001*step
variable p equal "(0.523*v_g+4.71225)*(57.29746936176986)"
variable t equal "sin(v_p)"
variable f equal "cos(v_p)"
fix zwalls1 all wall/gran/hertz/history primitive type 1 xplane 0.0
fix zwalls2 all wall/gran/hertz/history primitive type 1 xplane 2.4
fix cylwalls all wall/gran/hertz/history primitive type 1 xcylinder 0.157 0.157 0.157
region bc cylinder x 0.157 0.157 0.157 0.00 2.4 units box
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius constant 0.0125
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. 0. &
insert_every once overlapcheck yes all_in yes particles_in_region 3000 region bc
fix integr all nve/sphere
fix gravi all gravity 9.81 vector 0.0 v_f v_t
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
dump dmp all custom 800 dump.cilindro id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dmpp all atom 800 dump.gravity
#insert particles
run 100000
Actually it gives me a segmentation error, anyone who could help?
thank you for your attention
richti83 | Thu, 01/23/2014 - 10:18
fix insert
Hi,
I think fix insert/pack is not compatibel with variable gravity. Try to use an fixed gravity for settling for one step, after that run, unfix constant gravity and fix it in your way.
Let me know if this solves the problem.
Riccardo Maione | Thu, 01/23/2014 - 12:16
unluckily it doesn't work,
unluckily it doesn't work, here is my new script:
variable h equal step*0.000005
variable p equal "(0.523*h+4.71225)*(57.29746)"
variable m equal "sin(v_p)"
variable f equal "cos(v_f)"
lattice sc 0.1
atom_style granular
atom_modify map array
boundary p p p
echo screen
newton off
units si
communicate single vel yes
region reg block 0 2.4 0 0.314 0 0.314 units box
create_box 1 reg
neighbor 0.02 bin
neigh_modify delay 0
#proprietà dell'acciaio di prima simulazione
fix gravi all gravity 9.81 vector 0.0 -1.0 0.0
fix m1 all property/global youngsModulus peratomtype 2.1e11
fix m2 all property/global poissonsRatio peratomtype 0.27
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.95
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global coefficientRollingFriction peratomtypepair 1 0.3
#tipo di interazione, modello pseudo Hertz-Mindlin
pair_style gran/hertz/history rolling_friction cdt tangential_damping on
pair_coeff * *
#modello tipo Hertz-Mindlin, il timestep è di tipo 10^-6, visto che il tempo caratteristico di hertz è 10^-5
timestep 0.000005
fix zwalls1 all wall/gran/hertz/history primitive type 1 xplane 0.0
fix zwalls2 all wall/gran/hertz/history primitive type 1 xplane 2.4
fix cylwalls all wall/gran/hertz/history primitive type 1 xcylinder 0.157 0.157 0.157
#region of insertion
region bc cylinder x 0.157 0.157 0.157 0.00 2.4 units box
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius constant 0.0125
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.5 &
insert_every once overlapcheck yes all_in yes particles_in_region 3000 region bc
#apply nve integration to all particles
fix integr all nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
dump dmp all custom 1 dump.cilindro id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
#insert particles
run 100
unfix ins
unfix gravi
fix gravi all gravity 9.81 vector 0.0 v_f v_m
run 100
Riccardo Maione | Thu, 01/23/2014 - 15:07
it gives me the new error:
it gives me the new error: Invalid thermo keyword in variable formula (there is an error in the previous script the right one is "variable f equal "cos(v_p)" in place of variable f equal "cos(v_f)", I have already corrected it)
richti83 | Fri, 01/24/2014 - 18:45
1st: variable p equal (0.523
1st: variable p equal (0.523*v_h+4.71225)*(57.29746)
2nd: dont quote equal style variables, could be interpreted as text
3rd: I found a bug in fix_gravity, line 334: xvar,yvar,zvar are never evaluated, add these lines before double length=..
if (xstyle == EQUAL) xdir = input->variable->compute_equal(xvar);
if (ystyle == EQUAL) ydir = input->variable->compute_equal(yvar);
if (zstyle == EQUAL) zdir = input->variable->compute_equal(zvar);
My Inputscript look as this know and runs like a charm:
variable h equal step*0.000005
variable p equal (0.523*step*0.000005+4.71225)*(57.29746)
variable gx equal 0
variable gy equal sin(v_p)
variable gz equal cos(v_p)
lattice sc 0.1
atom_style granular
atom_modify map array
boundary p p p
echo screen
newton off
units si
communicate single vel yes
region reg block 0 2.4 0 0.314 0 0.314 units box
create_box 1 reg
neighbor 0.02 bin
neigh_modify delay 0
#proprietà dell'acciaio di prima simulazione
fix gravi all gravity 9.81 vector 0.0 -1.0 0.0
fix m1 all property/global youngsModulus peratomtype 2.1e11
fix m2 all property/global poissonsRatio peratomtype 0.27
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.95
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global coefficientRollingFriction peratomtypepair 1 0.3
#tipo di interazione, modello pseudo Hertz-Mindlin
pair_style gran/hertz/history rolling_friction cdt tangential_damping on
pair_coeff * *
#modello tipo Hertz-Mindlin, il timestep è di tipo 10^-6, visto che il tempo caratteristico di hertz è 10^-5
timestep 0.000005
fix zwalls1 all wall/gran/hertz/history primitive type 1 xplane 0.0
fix zwalls2 all wall/gran/hertz/history primitive type 1 xplane 2.4
fix cylwalls all wall/gran/hertz/history primitive type 1 xcylinder 0.157 0.157 0.157
#region of insertion
region bc cylinder x 0.157 0.157 0.157 0.00 2.4 units box
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 8000 radius constant 0.0125
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.5 &
insert_every once overlapcheck yes all_in yes particles_in_region 3000 region bc
#apply nve integration to all particles
fix integr all nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
dump dmp all custom 1 dump*.liggghts id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
#insert particles
run 100
unfix ins
unfix gravi
fix dbg all print 1 "${gx} ${gy} ${gz}" screen yes
fix gravi all gravity 9.81 vector v_gx v_gy v_gz
run 1000
Now you owe me a beer ;-)
Cheers,
Christian
techblog.richtisoft.de