simulation_shear.jl (8959B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import JLD 5 import PyPlot 6 import ArgParse 7 8 const render = false 9 10 function parse_command_line() 11 s = ArgParse.ArgParseSettings() 12 ArgParse.@add_arg_table s begin 13 "--E" 14 help = "Young's modulus [Pa]" 15 arg_type = Float64 16 default = 2e7 17 "--nu" 18 help = "Poisson's ratio [-]" 19 arg_type = Float64 20 default = 0.285 21 "--mu_d" 22 help = "dynamic friction coefficient [-]" 23 arg_type = Float64 24 default = 0.3 25 "--tensile_strength" 26 help = "the maximum tensile strength [Pa]" 27 arg_type = Float64 28 default = 0.0 29 "--r_min" 30 help = "minimum grain radius [m]" 31 arg_type = Float64 32 default = 5.0 33 "--r_max" 34 help = "maximum grain radius [m]" 35 arg_type = Float64 36 default = 1.35e3 37 default = 50.0 38 "--thickness" 39 help = "grain thickness [m]" 40 arg_type = Float64 41 default = 1.0 42 "--rotating" 43 help = "allow the grains to rotate" 44 arg_type = Bool 45 default = true 46 "--shearvel" 47 help = "shear velocity [m/s]" 48 arg_type = Float64 49 default = 1.0 50 "--seed" 51 help = "seed for the pseudo RNG" 52 arg_type = Int 53 default = 1 54 "id" 55 help = "simulation id" 56 required = true 57 end 58 return ArgParse.parse_args(s) 59 end 60 61 function report_args(parsed_args) 62 println("Parsed args:") 63 for (arg,val) in parsed_args 64 println(" $arg => $val") 65 end 66 end 67 68 function run_simulation(id::String, 69 E::Float64, 70 nu::Float64, 71 mu_d::Float64, 72 tensile_strength::Float64, 73 r_min::Float64, 74 r_max::Float64, 75 thickness::Float64, 76 rotating::Bool, 77 shearvel::Float64, 78 seed::Int) 79 80 ############################################################################ 81 #### SIMULATION PARAMETERS # 82 ############################################################################ 83 84 # Common simulation identifier 85 const id_prefix = id 86 87 # Shear velocity to apply to the top grains [m/s] 88 const vel_shear = shearvel 89 90 # Normal stresses for consolidation and shear [Pa] 91 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3] 92 93 # Simulation duration of individual steps [s] 94 const t_shear = 800.0 95 96 tau_p = Float64[] 97 tau_u = Float64[] 98 99 for N in N_list 100 101 ######################################################################## 102 #### Step 3: Shear the consolidated assemblage with a constant velocity# 103 ######################################################################## 104 105 sim_cons = Granular.createSimulation("$(id_prefix)-cons-N$(N)Pa") 106 sim = Granular.readSimulation(sim_cons) 107 108 sim.id = "$(id_prefix)-shear-N$(N)Pa-vel_shear$(vel_shear)m-s" 109 Granular.removeSimulationFiles(sim) 110 111 # Set all linear and rotational velocities to zero 112 Granular.zeroKinematics!(sim) 113 114 # Set current time to zero and reset output file counter 115 Granular.resetTime!(sim) 116 117 # Set the simulation time to run the shear experiment for 118 Granular.setTotalTime!(sim, t_shear) 119 120 y_bot = Inf 121 for grain in sim.grains 122 if y_bot > grain.lin_pos[2] - grain.contact_radius 123 y_bot = grain.lin_pos[2] - grain.contact_radius 124 end 125 Granular.disableOceanDrag!(grain) 126 end 127 # Run the shear experiment 128 time = Float64[] 129 shear_stress = Float64[] 130 shear_strain = Float64[] 131 dilation = Float64[] 132 const thickness_initial = sim.walls[1].pos - y_bot 133 x_min = +Inf 134 x_max = -Inf 135 for grain in sim.grains 136 if x_min > grain.lin_pos[1] - grain.contact_radius 137 x_min = grain.lin_pos[1] - grain.contact_radius 138 end 139 if x_max < grain.lin_pos[1] + grain.contact_radius 140 x_max = grain.lin_pos[1] + grain.contact_radius 141 end 142 end 143 const fixed_thickness = 2.0*r_max 144 const surface_area = (x_max - x_min) 145 shear_force = 0. 146 while sim.time < sim.time_total 147 148 # Prescribe the shear velocity to the uppermost grains 149 for grain in sim.grains 150 if grain.lin_pos[2] >= sim.walls[1].pos - fixed_thickness 151 grain.fixed = true 152 grain.color = 1 153 grain.allow_y_acc = true 154 grain.lin_vel[1] = vel_shear 155 elseif grain.lin_pos[2] > fixed_thickness # do not free bottom 156 grain.fixed = false 157 grain.color = 0 158 end 159 end 160 161 for i=1:100 162 Granular.run!(sim, single_step=true) 163 end 164 165 append!(time, sim.time) 166 167 # Determine the current shear stress 168 shear_force = 0.0 169 for grain in sim.grains 170 if grain.fixed && grain.allow_y_acc 171 shear_force += -grain.force[1] 172 end 173 end 174 append!(shear_stress, shear_force/surface_area) 175 176 # Determine the current shear strain 177 append!(shear_strain, sim.time*vel_shear/thickness_initial) 178 179 # Determine the current dilation 180 append!(dilation, (sim.walls[1].pos - y_bot)/thickness_initial) 181 182 end 183 184 # Try to render the simulation if `pvpython` is installed on the system 185 if render 186 Granular.render(sim, trim=false) 187 end 188 189 #Granular.writeSimulation(sim) 190 191 # Plot time vs. shear stress and dilation 192 PyPlot.subplot(211) 193 PyPlot.subplots_adjust(hspace=0.0) 194 ax1 = PyPlot.gca() 195 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 196 PyPlot.plot(time, shear_stress) 197 PyPlot.ylabel("Shear stress [Pa]") 198 PyPlot.subplot(212, sharex=ax1) 199 PyPlot.plot(time, dilation) 200 PyPlot.xlabel("Time [s]") 201 PyPlot.ylabel("Volumetric strain [-]") 202 PyPlot.savefig(sim.id * "-time_vs_shear-stress.pdf") 203 PyPlot.clf() 204 205 # Plot shear strain vs. shear stress and dilation 206 PyPlot.subplot(211) 207 PyPlot.subplots_adjust(hspace=0.0) 208 ax1 = PyPlot.gca() 209 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 210 PyPlot.plot(shear_strain, shear_stress) 211 PyPlot.ylabel("Shear stress [Pa]") 212 PyPlot.subplot(212, sharex=ax1) 213 PyPlot.plot(shear_strain, dilation) 214 PyPlot.xlabel("Shear strain [-]") 215 PyPlot.ylabel("Volumetric strain [-]") 216 PyPlot.savefig(sim.id * "-shear-strain_vs_shear-stress.pdf") 217 PyPlot.clf() 218 219 # Plot time vs. shear strain (boring when the shear experiment is rate 220 # controlled) 221 PyPlot.plot(time, shear_strain) 222 PyPlot.xlabel("Time [s]") 223 PyPlot.ylabel("Shear strain [-]") 224 PyPlot.savefig(sim.id * "-time_vs_shear-strain.pdf") 225 PyPlot.clf() 226 227 writedlm(sim.id * "-data.txt", [time, shear_strain, 228 shear_stress, dilation]) 229 230 append!(tau_p, maximum(shear_stress[1:100])) 231 append!(tau_u, mean(shear_stress[end-50:end])) 232 end 233 234 # Plot normal stress vs. peak/ultimate shear stress 235 writedlm(id_prefix * "-data.txt", [N_list, tau_p, tau_u]) 236 PyPlot.subplot(211) 237 PyPlot.subplots_adjust(hspace=0.0) 238 ax1 = PyPlot.gca() 239 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 240 PyPlot.plot(N_list, tau_p) 241 PyPlot.ylabel("Peak shear stress [Pa]") 242 PyPlot.subplot(212, sharex=ax1) 243 PyPlot.plot(N_list, tau_u) 244 PyPlot.xlabel("Shear strain [-]") 245 PyPlot.ylabel("Ultimate shear stress [Pa]") 246 PyPlot.xlim(0.-, maximum(N_list)*1.1) 247 PyPlot.savefig(id_prefix * "_mc.pdf") 248 PyPlot.clf() 249 end 250 251 function main() 252 parsed_args = parse_command_line() 253 report_args(parsed_args) 254 255 seed = parsed_args["seed"] 256 run_simulation(parsed_args["id"] * "-seed" * string(seed), 257 parsed_args["E"], 258 parsed_args["nu"], 259 parsed_args["mu_d"], 260 parsed_args["tensile_strength"], 261 parsed_args["r_min"], 262 parsed_args["r_max"], 263 parsed_args["thickness"], 264 parsed_args["rotating"], 265 parsed_args["shearvel"], 266 seed) 267 end 268 269 main()