###################################Header for General commands#######################################################
#particle packing by insertion and successive growing of particles
atom_style granular #Simulation of particles
atom_modify map array
boundary p f f #fixed boundaries -> particles will be deleted if leaving the simulation box and periodic boundaries in x axis
newton off #default
communicate single vel yes
units si
###################################System variables##################################################################
#Definition of boundaries
variable xmin_reg equal -1e-5
variable xmax_reg equal 3.0001e-3
variable ymin_reg equal -0.05001e-3
variable ymax_reg equal 8e-4
variable zmin_reg equal -1e-8
variable zmax_reg equal 3.0001e-3
variable xmin_inbox equal 0
variable xmax_inbox equal 3e-3
variable ymin_inbox equal 0
variable ymax_inbox equal 7.5e-4
variable zmin_inbox equal 0
variable zmax_inbox equal 3e-3
variable xmin_leftwall equal 1e-8
variable xmax_leftwall equal 2.999e-3
variable ymin_leftwall equal 7.499e-4
variable ymax_leftwall equal 7.99e-4
variable zmin_leftwall equal 1e-8
variable zmax_leftwall equal 2.99e-3
variable xmin_rightwall equal 1e-8
variable xmax_rightwall equal 2.999e-3
variable ymin_rightwall equal -0.049e-3
variable ymax_rightwall equal -1e-8
variable zmin_rightwall equal 1e-8
variable zmax_rightwall equal 2.999e-3
###################################Definition of simulationbox#######################################################
region reg block ${xmin_reg} ${xmax_reg} ${ymin_reg} ${ymax_reg} ${zmin_reg} ${zmax_reg} units box
region inbox block ${xmin_inbox} ${xmax_inbox} ${ymin_inbox} ${ymax_inbox} ${zmin_inbox} ${zmax_inbox} units box
region leftwall block ${xmin_leftwall} ${xmax_leftwall} ${ymin_leftwall} ${ymax_leftwall} ${zmin_leftwall} ${zmax_leftwall} units box
region rightwall block ${xmin_rightwall} ${xmax_rightwall} ${ymin_rightwall} ${ymax_rightwall} ${zmin_rightwall} ${zmax_rightwall} units box
create_box 1 reg
###################################neighbor list for element interaction######################################################
neighbor 0.0005 bin
neigh_modify delay 0 #default
###################################Definition of Material properties#################################################
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 7.e10
fix m2 all property/global poissonsRatio peratomtype 0.20
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global coefficientRollingFriction peratomtypepair 1 0
#################################Definition of the contact models####################################################
#New pair style
pair_style gran model hertz tangential history rolling_friction cdt #Hertzian without cohesion contact model
pair_coeff * * #default
###################################Generation and Loading of the Geometry .stl#######################################
fix ywall1 all wall/gran model hertz tangential history primitive type 1 yplane -0.05001e-3
fix ywall2 all wall/gran model hertz tangential history primitive type 1 yplane 8e-4
fix barrier1 all wall/gran model hertz tangential history primitive type 1 yplane 0
fix barrier2 all wall/gran model hertz tangential history primitive type 1 yplane 7.5e-4
fix xwall1 all wall/gran model hertz tangential history primitive type 1 xplane -1e-5
fix xwall2 all wall/gran model hertz tangential history primitive type 1 xplane 3.0001e-3
fix zwall all wall/gran model hertz tangential history primitive type 1 zplane -1e-8
fix zwall_out all wall/gran model hertz tangential history primitive type 1 zplane 3.0001e-3
#Definition of the timestep
timestep 1e-9
fix ts_check all check/timestep/gran 1000 0.2 0.2
###################################Generation and Insertion of the particles#########################################
#distributions for insertion
fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius constant 1e-4
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2650 radius constant 9e-5
fix pdd1 all particledistribution/discrete 32452843 2 pts1 1.0 pts2 0
fix pdd2 all particledistribution/discrete 32452867 2 pts1 0 pts2 1.0
###################################Specific variables for current simulation########################################
#parameters for gradually growing particle diameter
variable alphastart equal 0.40
variable alphatarget equal 0.65
variable growts equal 350000
variable growevery equal 40
variable relaxts equal 250000
#region and insertion
group nve_group region reg
fix ins nve_group insert/stream seed 49979687 distributiontemplate pdd1 vol ${alphastart} 200 region inbox
#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere
#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 100000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
print "Inserting Proppant Pack"
#insert the first particles
run 1
#dump dmp_paraview all custom 4000 ./post/dump_forces_limited_var_7.5e-4_1e03_*.liggghts id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
unfix ins
#undump dumpstl
#calculate grow rate
variable Rgrowrate equal (${alphatarget}/${alphastart})^(${growevery}/(3.*${growts}))
print "The radius grow rate is ${Rgrowrate}"
#do the diameter grow
compute rad all property/atom radius
variable d_grown atom ${Rgrowrate}*2.*c_rad
fix grow nve_group adapt ${growevery} atom diameter d_grown
neigh_modify every ${growevery} check no
neighbor 0.000005 bin
#run
run ${growts}
#let the packing relax
unfix grow
neigh_modify check yes
run ${relaxts}
variable pack equal ${growts}+${relaxts}
print "Creating Fracture Walls"
neighbor 0.0005 bin
lattice sc 5e-5
create_atoms 1 region leftwall units box
group leftwall_group region leftwall
set group leftwall_group diameter 5e-5 density 2650
create_atoms 1 region rightwall units box
group rightwall_group region rightwall
set group rightwall_group diameter 5e-5 density 2650
velocity rightwall_group zero linear
fix 2 rightwall_group nve/noforce
fix 3 rightwall_group freeze
#Granular Wall
fix 1 leftwall_group rigid group 1 leftwall_group force * on on off torque * off off off
neigh_modify exclude group leftwall_group leftwall_group
neigh_modify exclude group rightwall_group rightwall_group
neigh_modify exclude group leftwall_group rightwall_group
thermo 1000
#run 100000
unfix barrier
#save data
variable massleft equal mass(nve_group,inbox)
variable totalmass equal mass(nve_group,inbox)
variable averagevelocity equal vcm(nve_group,z,reg)
variable fracturewidth equal xcm(leftwall_group,y,reg)-xcm(rightwall_group,y,reg)
variable pressureleft equal fcm(leftwall_group,y,reg)
variable pressureright equal fcm(rightwall_group,y,reg)
fix screen all print 100000 "${massleft}, ${totalmass}, ${averagevelocity}, ${fracturewidth}, ${pressureleft}, ${pressureright}"
fix extra all print 4000 "${massleft}, ${totalmass}, ${averagevelocity}, ${fracturewidth}, ${pressureleft}, ${pressureright}" file flowingMAXs_203.txt
#relax
run 100000
#compact
print "Compacting Proppant Pack"
fix shmin leftwall_group addforce 0.0 0 0.0
run 400000
fix damp nve_group viscous 0.189
run 200000
unfix damp
unfix zwall_out
print "Applying Fracture Flow"
fix gr nve_group addforce 0. 0. 8.6167e-3
run 1000000
label loop
print "Entering mass loop"
variable masspres equal ${massleft}
run 40000
variable deltam equal v_masspres-v_massleft
if $deltam=<0.00000001 then "jump in.var_width_MAXp7.5e408.6167e3 continuing"
variable masspres delete
variable deltam delete
jump in.var_width_MAXp7.5e408.6167e3 loop
label continuing
run 100000
print "Done"