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 (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