Hi,
I am performing a biaxial compression simulation where I have a granular array surrounded by 4 extremely large granules to mimic a wall. In this simulation, I start off by fixing the left and bottom "walls" (big granules) and compress the array by adding a force (fix addforce) to move the top and right wall to compress the system. The force I use is to get a desired pressure in the confined state. Up to here all is good. I get a stress state which has the desired pressure. I add viscosity to dampen out the kinetic energy.
Next, I perform a test where I fix the top and right walls at the positions in the last step from the above explained simulation and set the viscosity to 0. I then run a NVE ensemble with all the walls fixed but in this compressed state. What I observe is the stress goes to 0 from the pressure I converge to in the confining simulation. I don't even have friction but I have no idea why a fully Hertzian system will decay to 0 stress after confining it. Is it because I am using a force control simulation? I don't know how i can get the desired stress state in a displacement control simulation by moving the walls but then again i notice that when I fix the walls, the stress dissipates to 0.
It would be greatly appreciated if you could let me know how to prevent this dissipation!
I have attached my input script below for reference. The pair style ep is pretty much hertzian contact for the settings I have provided! tw,bw,lw and rw stand for top wall, bottom wall, left and right wall. For reference, my wall (big granule) has a radius a billion times larger than the granules of my sample.
Many thanks!
DF
units si
atom_style sphere
atom_modify map array
boundary ff ff pp
newton off
comm_modify vel yes
dimension 2
log kn_84340800000.0_kt_57257000000.0_mu_0.5_ext_int3.lammps
neighbor 1e+6 bin
neigh_modify delay 0
neigh_modify exclude type 3 3
read_data data.sim
change_box all triclinic
variable kn equal 8.4248e+10
variable kn_wall equal 1.6849e+11
#pair_style hybrid gran/hertz/history ${kn} 57257000000.0 0 0 0.5 0 gran/ep/history 1.6849e+11 0 0 0 0 0 0.0 1 1 0 2
pair_style hybrid gran/hertz/history ${kn} 0 0 0 0.0 0 gran/ep/history 8.4248e+10 0 0 0 0 0 0.0 1 1 0 2
pair_coeff 1 1 gran/hertz/history
pair_coeff 1 2 gran/hertz/history
pair_coeff 2 2 gran/hertz/history
pair_coeff 1 3 gran/ep/history
pair_coeff 2 3 gran/ep/history
pair_coeff 3 3 gran/ep/history
group granules id 1:688
group walls id 689:692
group tw id 689
group bw id 690
group lw id 691
group rw id 692
group stationary id 690 691
group non_stationary id 1:688 689 692
group non_stationary_wall id 689 692
region atom_region block 0.04 0.36 0.04 1.96 -1e7 1e7
region atom_region2 block 0.12 0.28 0.8 1.8 -1e7 1e7
group atoms region atom_region
group atoms2 region atom_region2
timestep 1e-7
thermo 100
velocity stationary zero angular
velocity stationary zero linear
fix 1 all nve/sphere
fix 3 tw addforce 0 -16000 0
fix 6 rw addforce -80000 0 0
fix 7 granules viscous 500
fix 8 walls viscous 1000
fix 9 stationary freeze
fix 10 granules langevin 2 0 0.001 100000 omega yes
fix 11 tw setforce 0 NULL 0
fix 12 rw setforce NULL 0 0
compute peratom atoms stress/atom NULL pair
compute p atoms reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
compute p1 atoms reduce sum c_peratom[4] c_peratom[5] c_peratom[6]
variable press equal -(c_p[1]+c_p[2])/(2*0.04*(1.96-0.04)*(0.36-0.04))#((y[896]-1e+5)-(y[895]+1e+5))*((z[893]-1e+5)-(z[894]+1e+5)))
compute peratom2 atoms2 stress/atom NULL pair
compute p2_0 atoms2 reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]
compute p2_1 atoms2 reduce sum c_peratom2[4] c_peratom2[5] c_peratom2[6]
variable press2 equal -(c_p2_0[1]+c_p2_0[2])/(2*0.04*(0.28-0.12)*(1.8-0.8))
thermo_style custom step v_press ke v_press2 c_p[1] c_p[2] c_p[3] c_p1[1] c_p1[2] c_p1[3] c_p2_0[1] c_p2_0[2] c_p2_0[3] c_p2_1[1] c_p2_1[2] c_p2_1[3]
dump 3 all custom 20000 kn_84340800000.0_kt_57257000000.0_mu_0.5_ext_int3_no_friction_2.spheres id type mass radius x y vx vy omegaz ep es ef cj1 cj2 cj3 cj4 cj5 cj6 fx fy
dump_modify 3 format "%d %d %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e"
fix 2 all enforce2d
run 300000
unfix 7
unfix 3
unfix 6
unfix 8
unfix 1
fix 13 granules nve/sphere
fix 14 tw move linear 0 0 0
fix 15 rw move linear 0 0 0
fix 16 bw move linear 0 0 0
fix 17 lw move linear 0 0 0
run 100000
write_data kn_84340800000.0_kt_57257000000.0_mu_0.5_ext_int3_no_friction_2.sim