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