/
in.insert_particles
164 lines (133 loc) · 8.11 KB
/
in.insert_particles
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
## DEM Shear Test of wheat straw in LIGGGHTS
### Material Independent Properties
variable particle_diameter equal 2.50e-3 # particles diameter in meters
variable bond_out_diameter equal 2.50e-3 # fiber outer diameter in meters
variable bond_in_diameter equal 0.0 # fiber inner diameter in meters
variable bond_length equal 2.50e-3 # distance, in meters, between two particles in bond
variable bond_damp_type equal 1 # Use Yu Guo Damping
variable bond_damp_val equal 0.05 # NA
variable particle_density equal 2500.0 # kg/m3
variable fiber_contact_youngs_modulus equal 1.0e7 # Pa
variable wall_contact_youngs_modulus equal 180.0e9 # Pa
variable bond_youngs_modulus equal 1.0e7 # Pa
variable particle_poissons_ratio equal 0.3 # NA
variable wall_poissons_ratio equal 0.3 # NA
variable coef_res_pp equal 0.500 # particle-particle coefficient of restitution
variable coef_res_pw equal 0.300 # particle-wall coefficient of restitution
variable coef_res_ww equal 0.500 # wall-wall coefficient of restitution
variable coef_fri_pp equal 0.400 # particle-particle coefficient of friction
variable coef_fri_pw equal 0.600 # particle-wall coefficient of friction
variable coef_fri_ww equal 0.200 # wall-wall coefficient of friction
### Simulation Independent Parameters
variable ke_tol equal 1.0e-4 # Energy that we will run the simulation to obtain
variable res_tol equal 1.0e-6 # Change of height tolerance that must be met
variable shear_pressure equal 50000 # Pressure in N/m2
variable shear_box_length equal 0.101 # Length of shear box in meters
variable mass_shear_top_plate equal 0.401 # Mass of the top shear plate
variable restart_time equal 1.0e-2 # how often we write a restart file
variable fileprint_time equal 1.0e-2 # how often we print to the file in seconds
variable thermo_time equal 1.0e-3 # how often we print to the screen in seconds
variable output_time equal 1.0e-3 # how often data is written to csv file in seconds
### Material Dependent Properties
variable particle_radius equal 0.5*${particle_diameter}
variable bond_shear_modulus equal ${bond_youngs_modulus}/(2.0*(1.0+${particle_poissons_ratio}))
variable bond_out_per equal ${bond_out_diameter}/${particle_diameter}
variable bond_in_per equal ${bond_in_diameter}/${particle_diameter}
### Calculate dt using the bond model
variable r2 equal ${particle_radius}*${particle_radius}
variable r3 equal ${r2}*${particle_radius}
variable K equal ${bond_youngs_modulus}*PI*${r2}/${bond_length}
variable m equal 4.0*PI*${r3}*${particle_density}/3.0
variable w equal sqrt($K/$m)
variable dt equal 0.8/((1.0+3.0*${bond_damp_val})*$w)
### Simulation Dependent Parameters
variable restart_step equal ceil(${restart_time}/${dt})
variable fileprint_step equal ceil(${fileprint_time}/${dt})
variable thermo_step equal ceil(${thermo_time}/${dt})
variable output_step equal ceil(${output_time}/${dt})
variable shear_force equal ${shear_pressure}*(${shear_box_length}*${shear_box_length})+${mass_shear_top_plate}*9.81
## Build Simulation
# Set Atom Style
atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6
atom_modify map array
hard_particles yes
# Fix boundaries
boundary f f f
newton off
communicate single vel yes
units si
# Create simulation domain and tell liggghts that you have 2 types of particles
region domain block -0.125 0.225 -0.125 0.125 -0.01 0.125 units box
create_box 2 domain
# Use hertz-mindlin and granular fiber bonds
pair_style gran model hertz tangential history
bond_style gran
# Create a folder for restart files
shell mkdir restarts
# Set bin size and to check for a neighbor update every step
neighbor 0.001 bin
neigh_modify delay 0
# Set coefficients for contact and bond models
pair_coeff * * # none are needed
bond_coeff 1 ${bond_out_per} ${bond_in_per} ${bond_youngs_modulus} ${bond_shear_modulus} ${bond_damp_type} ${bond_damp_val} 1 1.0e32 1.0e32
# Set material properties
fix m1 all property/global youngsModulus peratomtype ${fiber_contact_youngs_modulus} ${wall_contact_youngs_modulus}
fix m2 all property/global poissonsRatio peratomtype ${particle_poissons_ratio} ${wall_poissons_ratio}
fix m3 all property/global coefficientRestitution peratomtypepair 2 ${coef_res_pp} ${coef_res_pw} &
${coef_res_pw} ${coef_res_ww}
fix m4 all property/global coefficientFriction peratomtypepair 2 ${coef_fri_pp} ${coef_fri_pw} &
${coef_fri_pw} ${coef_fri_ww}
# Import the geometry files and rotate/move to correct locations
fix BotShear all mesh/surface/stress file STL/ShearBox_TopPlate.STL type 2 &
scale 0.001 move 0.0 0.0 0.0 curvature_tolerant yes
# Top plate has a mass of 401 grams
fix TopShear all mesh/surface/stress/servo file STL/ShearBox_TopPlate.STL type 2 &
scale 0.001 rotate axis 1 0 0 angle 180 move 0.0 0.1015 0.090 &
curvature_tolerant yes com 0.05 0.05 0.08 ctrlPV force &
axis 0.0 0.0 -1.0 target_val ${shear_force} vel_max 5.0 &
kp 5.0e-2 ki 0.0 kd 0.0
fix BotHalfShear all mesh/surface/stress file STL/ShearBox_BottomHalf.STL &
type 2 scale 0.001 move 0.0 0.0 0.0 curvature_tolerant yes
fix TopHalfShear all mesh/surface/stress file STL/ShearBox_TopHalf.STL &
type 2 scale 0.001 move 0.0 0.0 0.025 curvature_tolerant yes
# Make all of the geometries able to interact with particles
fix wall all wall/gran model hertz tangential history mesh n_meshes 4 meshes &
BotShear TopShear BotHalfShear TopHalfShear
# Add primitive walls for particle insertion
fix wall_x1 all wall/gran model hertz tangential history primitive type 2 xplane 0.0
fix wall_x2 all wall/gran model hertz tangential history primitive type 2 xplane 0.1015
fix wall_y1 all wall/gran model hertz tangential history primitive type 2 yplane 0.0
fix wall_y2 all wall/gran model hertz tangential history primitive type 2 yplane 0.1015
# Add gravity
fix grav all gravity 9.81 vector 0.0 0.0 -1.0
# Integrate using velocity verlet
fix integr all nve/sphere
# Create our fibers
fix pts1 all particletemplate/multiplespheres 15485863 atom_type 1 &
density constant ${particle_density} nspheres 4 ntry 5000 spheres &
0.0000 0.0 0.0 0.00125 &
0.0025 0.0 0.0 0.00125 &
0.0050 0.0 0.0 0.00125 &
0.0075 0.0 0.0 0.00125 &
bonded no
# Set the distribution of particles
fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0
# Create a region to insert the particles
region fill_box block 0.0 0.1015 0.0 0.1015 0.0065 0.0835 units box
# Particle insertion
fix ins all insert/pack seed 32452867 distributiontemplate pdd1 &
maxattempt 50000 insert_every once overlapcheck yes orientation random &
all_in yes vel constant 0.0 0.0 -0.0 region fill_box &
particles_in_region 5000 ntry_mc 10000 check_dist_from_subdomain_border no
# Set the timestep to use
timestep ${dt}
# Set skin for for bond creation
variable bond_skin equal 1.000001*${particle_diameter}
# Bond the atoms together to form fibers
run 1
fix bondcr all bond/create/gran 1 1 1 ${bond_skin} 1 6 #every itype jtype cutoff btype newperts
run 1
fix_modify bondcr every 0 #do not create new bonds after this line
run 1
# Write the restart file
write_restart restarts/restart0.liggghts