ridging_simulation.jl (14405B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import PyPlot 5 import ArgParse 6 import DelimitedFiles 7 using Random 8 9 let 10 render = false 11 12 function parse_command_line() 13 s = ArgParse.ArgParseSettings() 14 ArgParse.@add_arg_table s begin 15 "--E" 16 help = "Young's modulus [Pa]" 17 arg_type = Float64 18 default = 2e7 19 "--nu" 20 help = "Poisson's ratio [-]" 21 arg_type = Float64 22 default = 0.285 23 "--mu_d" 24 help = "dynamic friction coefficient [-]" 25 arg_type = Float64 26 default = 0.3 27 "--tensile_strength" 28 help = "the maximum tensile strength [Pa]" 29 arg_type = Float64 30 default = 400e3 31 "--tensile_strength_std_dev" 32 help = "standard deviation of the maximum tensile strength [Pa]" 33 arg_type = Float64 34 default = 0.0 35 "--shear_strength" 36 help = "the maximum shear strength [Pa]" 37 arg_type = Float64 38 default = 200e3 39 "--r" 40 help = "grain radius [m]" 41 arg_type = Float64 42 default = 0.01 43 "--size_scaling" 44 help = "scale factor for r,nx1,nx2,ny1,ny2 [-]" 45 arg_type = Float64 46 default = 1.0 47 "--heal_rate" 48 help = "healing rate for new bonds [Pa/s]" 49 arg_type = Float64 50 default = 1.0 51 "--thickness" 52 help = "grain thickness [m]" 53 arg_type = Float64 54 default = 1. 55 "--rotating" 56 help = "allow the grains to rotate" 57 arg_type = Bool 58 default = true 59 "--compressive_velocity" 60 help = "compressive velocity [m/s]" 61 arg_type = Float64 62 default = 0.1 63 "--ny1" 64 help = "thickness in number of grains for left ice floe" 65 arg_type = Int 66 default = 10 67 "--ny2" 68 help = "thickness in number of grains for right ice floe" 69 arg_type = Int 70 default = 11 71 "--nx1" 72 help = "width in number of grains for left ice floe" 73 arg_type = Int 74 default = 100 75 "--nx2" 76 help = "width in number of grains for right ice floe" 77 arg_type = Int 78 default = 100 79 "--seed" 80 help = "seed for the pseudo RNG" 81 arg_type = Int 82 default = 1 83 "--one_floe" 84 help = "generate a single floe instead of two separate ones" 85 arg_type = Bool 86 default = false 87 "--many_floes" 88 help = "generate many floes instead of two separate ones" 89 arg_type = Bool 90 default = false 91 "id" 92 help = "simulation id" 93 required = true 94 end 95 return ArgParse.parse_args(s) 96 end 97 98 function report_args(parsed_args) 99 println("Parsed args:") 100 for (arg,val) in parsed_args 101 println(" $arg => $val") 102 end 103 end 104 105 function run_simulation(id::String, 106 E::Float64, 107 nu::Float64, 108 mu_d::Float64, 109 tensile_strength::Float64, 110 tensile_strength_std_dev::Float64, 111 shear_strength::Float64, 112 r::Float64, 113 size_scaling::Float64, 114 heal_rate::Float64, 115 thickness::Float64, 116 rotating::Bool, 117 compressive_velocity::Float64, 118 ny1::Int, 119 ny2::Int, 120 nx1::Int, 121 nx2::Int, 122 seed::Int, 123 one_floe::Bool=false, 124 many_floes::Bool=false) 125 126 ############################################################################ 127 #### SIMULATION PARAMETERS 128 ############################################################################ 129 130 # Common simulation identifier 131 id_prefix = id 132 133 # Grain mechanical properties 134 youngs_modulus = E # Elastic modulus [Pa] 135 poissons_ratio = nu # Shear-stiffness ratio [-] 136 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 137 138 # Simulation duration [s] 139 t_total = (nx1 + nx2)/2*r*2/compressive_velocity 140 141 # Temporal interval between output files 142 file_dt = t_total/200 143 144 # Misc parameters 145 water_density = 1000. 146 g = 9.8 147 #g = 0.0 148 y_water_surface = 0.0 149 150 r = r/size_scaling 151 nx1 = Int(round(nx1*size_scaling)) 152 nx2 = Int(round(nx2*size_scaling)) 153 ny1 = Int(round(ny1*size_scaling)) 154 ny2 = Int(round(ny2*size_scaling)) 155 156 Random.seed!(seed) 157 158 ############################################################################ 159 #### Create ice floes 160 ############################################################################ 161 sim = Granular.createSimulation("$(id_prefix)") 162 Granular.removeSimulationFiles(sim) 163 164 if one_floe 165 ny2 = ny1 166 Granular.regularPacking!(sim, [nx1+nx2, ny1], r, r, padding_factor=0., 167 tiling="triangular", 168 origo=[0.0, -0.66*ny1*r*2] 169 ) 170 # Add linear gradient in compressive velocity 171 max_x = (nx1+nx2)*r*2 172 for grain in sim.grains 173 grain.lin_vel[1] = (1.0 - grain.lin_pos[1]/max_x) * compressive_velocity 174 end 175 176 elseif many_floes 177 178 N_floes = 6 179 x = 0.0 180 for i=1:N_floes 181 182 t_total = nx1*N_floes*r/compressive_velocity 183 184 nx = nx1 + Int(round(nx1*0.5*rand())) 185 ny = ny1 + Int(round(ny1*0.5*rand())) 186 187 # Left ice floe 188 Granular.regularPacking!(sim, [nx, ny], r, r, padding_factor=0., 189 tiling="triangular", 190 origo=[x, -0.66*ny1*r*2] 191 ) 192 193 if i == 1 # impose velocity onty first ice floe 194 for grain in sim.grains 195 grain.lin_vel[1] = compressive_velocity 196 end 197 end 198 199 x += nx * 2.0*r + 4.0*r 200 end 201 202 else # two ice floes 203 # Left ice floe 204 Granular.regularPacking!(sim, [nx1, ny1], r, r, padding_factor=0., 205 tiling="triangular", 206 origo=[0.0, -0.66*ny1*r*2] 207 ) 208 for grain in sim.grains 209 grain.lin_vel[1] = compressive_velocity 210 end 211 212 # Right ice floe 213 Granular.regularPacking!(sim, [nx2, ny2], r, r, padding_factor=0., 214 tiling="triangular", 215 origo=[nx1*2*r + 10*r*size_scaling, 216 -0.66*ny2*r*2] 217 ) 218 end 219 220 for grain in sim.grains 221 grain.contact_radius *= 1.0 + 1e-6 222 grain.areal_radius *= 1.0 + 1e-6 223 end 224 225 if many_floes 226 Granular.fitGridToGrains!(sim, sim.ocean, padding=r*100, verbose=true) 227 228 else 229 # Create a grid for contact searching spanning the extent of the grains 230 sim.ocean = Granular.createRegularOceanGrid([(nx1+nx2) - 10, 231 max(ny1,ny2)*10 - 2, 1], 232 [(nx1+nx2)*r*2 + 11*r, 233 max(ny1,ny2)*2*r*10, 234 2*r], 235 origo=[0., -max(ny1,ny2)*2*r*8*0.75]) 236 end 237 238 # Make the top and bottom boundaries impermeable, and the side boundaries 239 # periodic, which will come in handy during shear 240 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable") 241 242 # Set grain mechanical properties 243 for grain in sim.grains 244 grain.youngs_modulus = youngs_modulus 245 grain.poissons_ratio = poissons_ratio 246 grain.tensile_strength = abs(tensile_strength + randn()*tensile_strength_std_dev) 247 grain.shear_strength = abs(shear_strength + randn()*tensile_strength_std_dev) 248 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 249 grain.contact_dynamic_friction = contact_dynamic_friction 250 grain.rotating = rotating 251 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 252 end 253 254 Granular.setTimeStep!(sim) 255 Granular.setTotalTime!(sim, t_total) 256 Granular.setOutputFileInterval!(sim, file_dt) 257 if render 258 Granular.plotGrains(sim, verbose=true) 259 end 260 261 grainArrays = Granular.convertGrainDataToArrays(sim) 262 x_max = maximum(grainArrays.lin_pos[1,:]) 263 264 # fix velocities of outermost parts of the ice floes 265 if many_floes 266 fix_dist = nx1*r*2 267 else 268 fix_dist = r*10 269 end 270 271 for grain in sim.grains 272 if grain.lin_pos[1] < fix_dist 273 grain.fixed = true 274 grain.allow_y_acc = true 275 276 elseif grain.lin_pos[1] > x_max - fix_dist 277 grain.fixed = true 278 grain.allow_y_acc = true 279 end 280 end 281 282 ############################################################################ 283 #### RUN THE SIMULATION 284 ############################################################################ 285 theta = 0.0; submerged_volume = 0.0 286 time = Float64[] 287 N = Float64[] # normal stress on leftmost fixed particles 288 τ = Float64[] # shear stress on leftmost fixed particles 289 A = sim.grains[1].thickness * ny1*r*2 290 while sim.time < sim.time_total 291 292 append!(time, sim.time) 293 294 # Record cumulative force on leftmost fixed particles 295 normal_force = 0.0 296 tangential_force = 0.0 297 for grain in sim.grains 298 if grain.fixed && grain.lin_pos[1] < 0.7*x_max 299 normal_force += grain.force[1] 300 tangential_force += grain.force[2] 301 end 302 end 303 append!(N, -normal_force/A) # compressive = positive 304 append!(τ, tangential_force/A) 305 306 # Run for 100 time steps before measuring bulk forces 307 for i=1:100 308 309 if sim.time_iteration < 3 310 # Age (and strengthen) all existing bonds 311 for grain in sim.grains 312 for ic=1:length(grain.contact_age) 313 grain.contact_age[ic] = 1e16 314 end 315 end 316 elseif sim.time_iteration == 3 317 # weaken any new bonding 318 for grain in sim.grains 319 grain.strength_heal_rate = heal_rate # new bond stengthening 320 end 321 end 322 323 # Adjust body forces, assuming water surface at y=0 324 for grain in sim.grains 325 326 # Gravitational force 327 Granular.setBodyForce!(grain, [0.0, -grain.mass*g]) 328 329 330 if grain.lin_pos[2] + grain.areal_radius < y_water_surface 331 # Add buoyant force, fully submerged 332 Granular.addBodyForce!(grain, [0.0, 333 water_density*grain.volume*g]) 334 # set ocean drag coefficient 335 grain.ocean_drag_coeff_vert = 0.85 336 337 338 elseif grain.lin_pos[2] - grain.areal_radius < y_water_surface 339 # Add buoyant force, partially submerged 340 341 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_radius) 342 if grain.lin_pos[2] < y_water_surface 343 submerged_volume = grain.volume - 344 (grain.areal_radius^2/2* 345 (theta - sin(theta))*grain.thickness) 346 else 347 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_radius) 348 submerged_volume = grain.areal_radius^2/2* 349 (theta - sin(theta))*grain.thickness 350 end 351 352 Granular.addBodyForce!(grain, 353 [0.0, water_density*submerged_volume*g]) 354 grain.ocean_drag_coeff_vert = 0.85 355 else 356 # above water 357 grain.ocean_drag_coeff_vert = 0.0 358 end 359 grain.color = Int(round(grain.external_body_force[2]*1000)) 360 end 361 Granular.run!(sim, single_step=true) 362 end 363 364 end 365 366 # Plot normal stress and shear stress vs time 367 DelimitedFiles.writedlm(id_prefix * "-data.txt", [time, N, τ]) 368 PyPlot.subplot(211) 369 PyPlot.subplots_adjust(hspace=0.0) 370 ax1 = PyPlot.gca() 371 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 372 PyPlot.plot(time, N/1e3) 373 PyPlot.ylabel("Compressive stress [kPa]") 374 PyPlot.subplot(212, sharex=ax1) 375 PyPlot.plot(time, τ/1e3) 376 PyPlot.xlabel("Time [s]") 377 PyPlot.ylabel("Shear stress [kPa]") 378 PyPlot.tight_layout() 379 PyPlot.savefig(sim.id * "-time_vs_stress.pdf") 380 PyPlot.clf() 381 382 # Create GSD plot to signify that simulation is complete 383 #Granular.writeSimulation(sim) 384 #if render 385 # Granular.plotGrainSizeDistribution(sim) 386 #end 387 end 388 389 function main() 390 parsed_args = parse_command_line() 391 report_args(parsed_args) 392 393 seed = parsed_args["seed"] 394 run_simulation(parsed_args["id"] * "-seed" * string(seed), 395 parsed_args["E"], 396 parsed_args["nu"], 397 parsed_args["mu_d"], 398 parsed_args["tensile_strength"], 399 parsed_args["tensile_strength_std_dev"], 400 parsed_args["shear_strength"], 401 parsed_args["r"], 402 parsed_args["size_scaling"], 403 parsed_args["heal_rate"], 404 parsed_args["thickness"], 405 parsed_args["rotating"], 406 parsed_args["compressive_velocity"], 407 parsed_args["ny1"], 408 parsed_args["ny2"], 409 parsed_args["nx1"], 410 parsed_args["nx2"], 411 seed, 412 parsed_args["one_floe"], 413 parsed_args["many_floes"]) 414 end 415 416 main() 417 end