tridging_simulation.jl - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tridging_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