Inconsistencies during uniaxial compression

Submitted by Dario Passafiume on Tue, 07/27/2021 - 13:04

Hi everyone,

I am trying to perform a quasi-static uniaxial compression of an assembly made of a cubic box with inside 5000 particles randomly distributed with a radius of 0.00025 meters.

I am realizing this by using the "fix deform" command and a double loop. The external loop is for compression, the internal one for dissipation of kinetic energy realized with the "fix viscous" command.

I want to get stress-strain curves for the whole assembly.
To compute the overall stress tensor I am using the "compute pressure" command. However, I am getting a strange behavior at the beginning of the simulation; in particular, the mean kinetic energy (so the particles' velocities as well) and the stresses (so the particles' forces as well) remain equal to zero up to a strain of -0.44%. This is of course unphysical because particles are already moving and compressing against each other, as I verified in the post-processing with Paraview.

I tried to perform a simpler calculation without the use of loops, so just compressing with "fix deform" but the behavior is exactly the same up to that value of the strain. So I think there might be something wrong either with the way I am evaluating the stresses or with my geometry or even the way I'm applying the deformation.

I thought of evaluating the stress tensor with the "fix ave/euler" command but I think this is not helpful since it calculates the tensor in cells the system is divided in. Instead what I need is the tensor of my whole assembly.

Can you help me in figuring out what might be the reason for that behavior, please?

Here below you find my input script.

####UNIAXIAL COMPRESSION TEST#####

atom_style granular
boundary p p p
newton off
echo both
units si

hard_particles yes

communicate single vel yes

read_data KIT_DEM_CODE_data

neighbor 0.0001 bin
neigh_modify delay 0 page 500000 one 20000

atom_modify sort 0 0

#####Material properties

fix m1 all property/global youngsModulus peratomtype 90.e9
fix m2 all property/global poissonsRatio peratomtype 0.25
fix m3 all property/global coefficientRestitution peratomtypepair 1 1
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1

#####New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion and rolling
pair_coeff * *

timestep 0.000000001
variable dt equal step

#timestep check
#fix ts_check all check/timestep/gran 5 0.2 0.2 warn yes

####integrator
fix 1 all nve/sphere

####compute overlap between particles
compute 1 all pair/gran/local id delta
#save overlap of particles throughout the simulation
dump overlap all local 10 overlap/overlap*.txt c_1[4]

#####compute kinetic energy of the system (summation of all particles) and the average value
compute 2 all ke update_on_run_end yes
compute 3 all erotate update_on_run_end yes
variable KE1 equal c_2
variable KE2 equal c_3
#variable AVE_kin_en equal 1e-9
variable AVEKE equal (v_KE1+v_KE2)/5000

####compute overall stress tensor of the assembly and z component for the stress-strain curve
compute Temp all temp
compute 4 all pressure Temp virial
variable stressXX equal c_4[1]
variable stressYY equal c_4[2]
variable stressZZ equal c_4[3]
variable trace equal abs(c_4[1]+c_4[2]+c_4[3])

#### Store final cell length for strain calculations
variable tmp equal "lz"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"

####output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms v_AVEKE v_stressZZ
thermo 1
thermo_modify lost ignore norm no

run 0
dump dmp all custom/vtk 1 post/packing_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
fix 5 all print 1 "${dt} $s ${stressZZ}" file Stress_Strain.txt

####Applying uniaxial compression (loading cicle)
label loopa
variable a loop 500
fix press all deform 1 z delta 0 -2.61820e-7 # negative strain application along z axis
variable strain equal (lz-${L0})/${L0}
variable s equal v_strain
run 1
label loopb
variable b loop 100
fix press all deform 1 z delta 0 -2.61820e-7
unfix press
fix 3 all viscous 20
if "${AVEKE} < 1e-12" then "jump in.packing_Xia_July19 break1"
run 1
next b
jump in.packing_Xia_July19 loopb
label break1
variable b delete
unfix 3
next a
jump in.packing_Xia_July19 loopa

print "ALL DONE"