 |
tridging_simulation.jl - seaice-experiments - sea ice experiments using Granula… |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
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 [… |
|
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_fact… |
|
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) * compress… |
|
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… |
|
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… |
|
271 |
|
272 else |
|
273 # Create a grid for contact searching spanning the extent of the… |
|
274 #sim.ocean = Granular.createRegularOceanGrid([(nx1+nx2) - 10, |
|
275 # max(ny1,ny2)*10 - … |
|
276 # [(nx1+nx2)*r*2 + 11… |
|
277 # max(ny1,ny2)*2*r*1… |
|
278 # 2*r], |
|
279 # origo=[0., -max(ny1… |
|
280 Granular.fitGridToGrains!(sim, sim.ocean, padding=r*100, verbose… |
|
281 end |
|
282 |
|
283 # Make the top and bottom boundaries impermeable, and the side bound… |
|
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_… |
|
292 grain.shear_strength = abs(shear_strength + randn()*tensile_stre… |
|
293 grain.contact_static_friction = contact_dynamic_friction # mu_s… |
|
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 sten… |
|
362 end |
|
363 end |
|
364 |
|
365 # remove velocity constraint on grains that have left the fi… |
|
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 … |
|
372 grain.fixed = true |
|
373 grain.allow_y_acc = true |
|
374 else |
|
375 if grain.fixed == true |
|
376 newgrain = deepc… |
|
377 grain.fixed = fa… |
|
378 # keep shifting … |
|
379 if grain.lin_vel… |
|
380 while !G… |
|
381 … |
|
382 … |
|
383 … |
|
384 end |
|
385 else |
|
386 while !G… |
|
387 … |
|
388 … |
|
389 … |
|
390 end |
|
391 end |
|
392 newgrain.lin_dis… |
|
393 newgrain.n_conta… |
|
394 newgrain.contact… |
|
395 Granular.addGrai… |
|
396 i_newgrain = len… |
|
397 println("freeing… |
|
398 println("adding … |
|
399 Granular.sortGra… |
|
400 Granular.findCon… |
|
401 println("new gra… |
|
402 # max out contac… |
|
403 for g2 in sim.gr… |
|
404 if i_new… |
|
405 … |
|
406 … |
|
407 … |
|
408 … |
|
409 … |
|
410 … |
|
411 … |
|
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_surfa… |
|
427 # Add buoyant force, fully submerged |
|
428 Granular.addBodyForce!(grain, [0.0, |
|
429 water_density*grain.v… |
|
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_s… |
|
435 # Add buoyant force, partially submerged |
|
436 |
|
437 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_rad… |
|
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… |
|
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… |
|
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 |
