ridging_simulation.jl (17625B)
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 "--continuous_ice" 92 help = "add ice continuously and compress from both sides" 93 arg_type = Bool 94 default = false 95 "id" 96 help = "simulation id" 97 required = true 98 end 99 return ArgParse.parse_args(s) 100 end 101 102 function report_args(parsed_args) 103 println("Parsed args:") 104 for (arg,val) in parsed_args 105 println(" $arg => $val") 106 end 107 end 108 109 function run_simulation(id::String, 110 E::Float64, 111 nu::Float64, 112 mu_d::Float64, 113 tensile_strength::Float64, 114 tensile_strength_std_dev::Float64, 115 shear_strength::Float64, 116 r::Float64, 117 size_scaling::Float64, 118 heal_rate::Float64, 119 thickness::Float64, 120 rotating::Bool, 121 compressive_velocity::Float64, 122 ny1::Int, 123 ny2::Int, 124 nx1::Int, 125 nx2::Int, 126 seed::Int, 127 one_floe::Bool=false, 128 many_floes::Bool=false, 129 continuous_ice::Bool=false) 130 131 ############################################################################ 132 #### SIMULATION PARAMETERS 133 ############################################################################ 134 135 # Common simulation identifier 136 id_prefix = id 137 138 # Grain mechanical properties 139 youngs_modulus = E # Elastic modulus [Pa] 140 poissons_ratio = nu # Shear-stiffness ratio [-] 141 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 142 143 # Simulation duration [s] 144 t_total = (nx1 + nx2)/2*r*2/compressive_velocity 145 146 # Temporal interval between output files 147 file_dt = t_total/200 148 149 # Misc parameters 150 water_density = 1000. 151 g = 9.8 152 #g = 0.0 153 y_water_surface = 0.0 154 155 r = r/size_scaling 156 nx1 = Int(round(nx1*size_scaling)) 157 nx2 = Int(round(nx2*size_scaling)) 158 ny1 = Int(round(ny1*size_scaling)) 159 ny2 = Int(round(ny2*size_scaling)) 160 161 Random.seed!(seed) 162 163 ############################################################################ 164 #### Create ice floes 165 ############################################################################ 166 sim = Granular.createSimulation("$(id_prefix)") 167 Granular.removeSimulationFiles(sim) 168 169 if one_floe 170 ny2 = ny1 171 Granular.regularPacking!(sim, [nx1+nx2, ny1], r, r, padding_factor=0., 172 tiling="triangular", 173 origo=[0.0, -0.66*ny1*r*2] 174 ) 175 # Add linear gradient in compressive velocity 176 max_x = (nx1+nx2)*r*2 177 for grain in sim.grains 178 grain.lin_vel[1] = (1.0 - grain.lin_pos[1]/max_x) * compressive_velocity 179 end 180 181 elseif many_floes 182 183 N_floes = 6 184 x = 0.0 185 for i=1:N_floes 186 187 t_total = nx1*N_floes*r/compressive_velocity 188 189 nx = nx1 + Int(round(nx1*0.5*rand())) 190 ny = ny1 + Int(round(ny1*0.5*rand())) 191 192 N0 = length(sim.grains) 193 # Left ice floe 194 Granular.regularPacking!(sim, [nx, ny], r, r, padding_factor=0., 195 tiling="triangular", 196 origo=[x, -0.66*ny1*r*2] 197 ) 198 if i == 1 199 for grain in sim.grains 200 grain.lin_vel[1] = 0.5*compressive_velocity 201 end 202 end 203 204 N1 = length(sim.grains) 205 206 for in=(N0 + 1):N1 207 sim.grains[in].color = i 208 end 209 210 if i == N_floes 211 for in=(N0 + 1):N1 212 sim.grains[in].lin_vel[1] = -0.5*compressive_velocity 213 end 214 end 215 216 x += nx * 2.0*r + 4.0*r 217 end 218 219 elseif continuous_ice 220 221 # Left ice floe 222 Granular.regularPacking!(sim, [nx1, ny1], r, r, padding_factor=0., 223 tiling="triangular", 224 origo=[0.0, -0.66*ny1*r*2] 225 ) 226 for grain in sim.grains 227 grain.lin_vel[1] = 0.5*compressive_velocity 228 end 229 230 # Right ice floe 231 Granular.regularPacking!(sim, [nx2, ny2], r, r, padding_factor=0., 232 tiling="triangular", 233 origo=[nx1*2*r + 10*r*size_scaling, 234 -0.66*ny2*r*2] 235 ) 236 for grain in sim.grains 237 if grain.lin_vel[1] == 0.0 238 grain.lin_vel[1] = -0.5*compressive_velocity 239 end 240 end 241 242 else # two ice floes 243 244 # Left ice floe 245 Granular.regularPacking!(sim, [nx1, ny1], r, r, padding_factor=0., 246 tiling="triangular", 247 origo=[0.0, -0.66*ny1*r*2] 248 ) 249 grainArrays = Granular.convertGrainDataToArrays(sim) 250 x_min = minimum(grainArrays.lin_pos[1,:]) 251 for grain in sim.grains 252 grain.lin_vel[1] = compressive_velocity 253 grain.lin_pos[1] -= x_min 254 end 255 256 # Right ice floe 257 Granular.regularPacking!(sim, [nx2, ny2], r, r, padding_factor=0., 258 tiling="triangular", 259 origo=[nx1*2*r + 10*r*size_scaling, 260 -0.66*ny2*r*2] 261 ) 262 end 263 264 for grain in sim.grains 265 grain.contact_radius *= 1.0 + 1e-6 266 grain.areal_radius *= 1.0 + 1e-6 267 end 268 269 if many_floes 270 Granular.fitGridToGrains!(sim, sim.ocean, padding=r*100, verbose=true) 271 272 else 273 # Create a grid for contact searching spanning the extent of the grains 274 #sim.ocean = Granular.createRegularOceanGrid([(nx1+nx2) - 10, 275 # max(ny1,ny2)*10 - 2, 1], 276 # [(nx1+nx2)*r*2 + 11*r, 277 # max(ny1,ny2)*2*r*10, 278 # 2*r], 279 # origo=[0., -max(ny1,ny2)*2*r*8*0.75]) 280 Granular.fitGridToGrains!(sim, sim.ocean, padding=r*100, verbose=true) 281 end 282 283 # Make the top and bottom boundaries impermeable, and the side boundaries 284 # periodic, which will come in handy during shear 285 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable") 286 287 # Set grain mechanical properties 288 for grain in sim.grains 289 grain.youngs_modulus = youngs_modulus 290 grain.poissons_ratio = poissons_ratio 291 grain.tensile_strength = abs(tensile_strength + randn()*tensile_strength_std_dev) 292 grain.shear_strength = abs(shear_strength + randn()*tensile_strength_std_dev) 293 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 294 grain.contact_dynamic_friction = contact_dynamic_friction 295 grain.rotating = rotating 296 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 297 end 298 299 Granular.setTimeStep!(sim) 300 Granular.setTotalTime!(sim, t_total) 301 Granular.setOutputFileInterval!(sim, file_dt) 302 if render 303 Granular.plotGrains(sim, verbose=true) 304 end 305 306 grainArrays = Granular.convertGrainDataToArrays(sim) 307 x_max = maximum(grainArrays.lin_pos[1,:]) 308 x_min = minimum(grainArrays.lin_pos[1,:]) 309 310 # fix velocities of outermost parts of the ice floes 311 if many_floes 312 fix_dist = nx1*r*2 313 else 314 fix_dist = r*10 315 end 316 317 for grain in sim.grains 318 if abs(grain.lin_vel[1]) > 0.0 319 grain.fixed = true 320 grain.allow_y_acc = true 321 end 322 end 323 324 ############################################################################ 325 #### RUN THE SIMULATION 326 ############################################################################ 327 theta = 0.0; submerged_volume = 0.0 328 time = Float64[] 329 N = Float64[] # normal stress on leftmost fixed particles 330 τ = Float64[] # shear stress on leftmost fixed particles 331 A = sim.grains[1].thickness * ny1*r*2 332 while sim.time < sim.time_total 333 334 append!(time, sim.time) 335 336 # Record cumulative force on leftmost fixed particles 337 normal_force = 0.0 338 tangential_force = 0.0 339 for grain in sim.grains 340 if grain.fixed && grain.lin_pos[1] < 0.7*x_max 341 normal_force += grain.force[1] 342 tangential_force += grain.force[2] 343 end 344 end 345 append!(N, -normal_force/A) # compressive = positive 346 append!(τ, tangential_force/A) 347 348 # Run for 100 time steps before measuring bulk forces 349 for i=1:100 350 351 if sim.time_iteration < 3 352 # Age (and strengthen) all existing bonds 353 for grain in sim.grains 354 for ic=1:length(grain.contact_age) 355 grain.contact_age[ic] = 1e16 356 end 357 end 358 elseif sim.time_iteration == 3 359 # weaken any new bonding 360 for grain in sim.grains 361 grain.strength_heal_rate = heal_rate # new bond stengthening 362 end 363 end 364 365 # remove velocity constraint on grains that have left the fixed zone 366 if continuous_ice 367 for grain in sim.grains 368 if grain.lin_pos[1] < fix_dist 369 grain.fixed = true 370 grain.allow_y_acc = true 371 elseif grain.lin_pos[1] > x_max - fix_dist 372 grain.fixed = true 373 grain.allow_y_acc = true 374 else 375 if grain.fixed == true 376 newgrain = deepcopy(grain) 377 grain.fixed = false 378 # keep shifting outwards until free space is found 379 if grain.lin_vel[1] > 0.0 380 while !Granular.checkForContacts(sim, sim.ocean, 381 newgrain.lin_pos, r*0.1, 382 return_when_overlap_found=true) 383 newgrain.lin_pos[1] -= 2.0*r 384 end 385 else 386 while !Granular.checkForContacts(sim, sim.ocean, 387 newgrain.lin_pos, r*0.1, 388 return_when_overlap_found=true) 389 newgrain.lin_pos[1] += 2.0*r 390 end 391 end 392 newgrain.lin_disp .= zeros(3) 393 newgrain.n_contacts = 0 394 newgrain.contacts .= 0 395 Granular.addGrain!(sim, newgrain) 396 i_newgrain = length(sim.grains) 397 println("freeing grain at $(grain.lin_pos)") 398 println("adding new grain $(i_newgrain) at $(newgrain.lin_pos)") 399 Granular.sortGrainsInGrid!(sim, sim.ocean, verbose=false) 400 Granular.findContactsInGrid!(sim, sim.ocean) 401 println("new grain n_contacts: $(sim.grains[end].n_contacts)") 402 # max out contact age in list for opposite grain 403 for g2 in sim.grains 404 if i_newgrain in g2.contacts 405 println("contact list: $(g2.contacts)") 406 for ic=1:g2.n_contacts 407 if g2.contacts[ic] == i_newgrain 408 g2.contact_age[ic] = 1e16 409 end 410 end 411 println("contact ages: $(g2.contact_age)") 412 end 413 end 414 end 415 end 416 end 417 end 418 419 # Adjust body forces, assuming water surface at y=0 420 for grain in sim.grains 421 422 # Gravitational force 423 Granular.setBodyForce!(grain, [0.0, -grain.mass*g]) 424 425 426 if grain.lin_pos[2] + grain.areal_radius < y_water_surface 427 # Add buoyant force, fully submerged 428 Granular.addBodyForce!(grain, [0.0, 429 water_density*grain.volume*g]) 430 # set ocean drag coefficient 431 grain.ocean_drag_coeff_vert = 0.85 432 433 434 elseif grain.lin_pos[2] - grain.areal_radius < y_water_surface 435 # Add buoyant force, partially submerged 436 437 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_radius) 438 if grain.lin_pos[2] < y_water_surface 439 submerged_volume = grain.volume - 440 (grain.areal_radius^2/2* 441 (theta - sin(theta))*grain.thickness) 442 else 443 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_radius) 444 submerged_volume = grain.areal_radius^2/2* 445 (theta - sin(theta))*grain.thickness 446 end 447 448 Granular.addBodyForce!(grain, 449 [0.0, water_density*submerged_volume*g]) 450 grain.ocean_drag_coeff_vert = 0.85 451 else 452 # above water 453 grain.ocean_drag_coeff_vert = 0.0 454 end 455 end 456 Granular.run!(sim, single_step=true) 457 end 458 459 end 460 461 # Plot normal stress and shear stress vs time 462 DelimitedFiles.writedlm(id_prefix * "-data.txt", [time, N, τ]) 463 PyPlot.subplot(211) 464 PyPlot.subplots_adjust(hspace=0.0) 465 ax1 = PyPlot.gca() 466 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 467 PyPlot.plot(time, N/1e3) 468 PyPlot.ylabel("Compressive stress [kPa]") 469 PyPlot.subplot(212, sharex=ax1) 470 PyPlot.plot(time, τ/1e3) 471 PyPlot.xlabel("Time [s]") 472 PyPlot.ylabel("Shear stress [kPa]") 473 PyPlot.tight_layout() 474 PyPlot.savefig(sim.id * "-time_vs_stress.pdf") 475 PyPlot.clf() 476 477 # Create GSD plot to signify that simulation is complete 478 #Granular.writeSimulation(sim) 479 #if render 480 # Granular.plotGrainSizeDistribution(sim) 481 #end 482 end 483 484 function main() 485 parsed_args = parse_command_line() 486 report_args(parsed_args) 487 488 seed = parsed_args["seed"] 489 run_simulation(parsed_args["id"] * "-seed" * string(seed), 490 parsed_args["E"], 491 parsed_args["nu"], 492 parsed_args["mu_d"], 493 parsed_args["tensile_strength"], 494 parsed_args["tensile_strength_std_dev"], 495 parsed_args["shear_strength"], 496 parsed_args["r"], 497 parsed_args["size_scaling"], 498 parsed_args["heal_rate"], 499 parsed_args["thickness"], 500 parsed_args["rotating"], 501 parsed_args["compressive_velocity"], 502 parsed_args["ny1"], 503 parsed_args["ny2"], 504 parsed_args["nx1"], 505 parsed_args["nx2"], 506 seed, 507 parsed_args["one_floe"], 508 parsed_args["many_floes"], 509 parsed_args["continuous_ice"]) 510 end 511 512 main() 513 end