 |
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 (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 [… |
|
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_fact… |
|
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) * compress… |
|
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… |
|
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… |
|
227 |
|
228 else |
|
229 # Create a grid for contact searching spanning the extent of the… |
|
230 sim.ocean = Granular.createRegularOceanGrid([(nx1+nx2) - 10, |
|
231 max(ny1,ny2)*10 - 2… |
|
232 [(nx1+nx2)*r*2 + 11*… |
|
233 max(ny1,ny2)*2*r*10, |
|
234 2*r], |
|
235 origo=[0., -max(ny1,… |
|
236 end |
|
237 |
|
238 # Make the top and bottom boundaries impermeable, and the side bound… |
|
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_… |
|
247 grain.shear_strength = abs(shear_strength + randn()*tensile_stre… |
|
248 grain.contact_static_friction = contact_dynamic_friction # mu_s… |
|
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 sten… |
|
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_surfa… |
|
331 # Add buoyant force, fully submerged |
|
332 Granular.addBodyForce!(grain, [0.0, |
|
333 water_density*grain.v… |
|
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_s… |
|
339 # Add buoyant force, partially submerged |
|
340 |
|
341 theta = 2*acos(abs(grain.lin_pos[2])/grain.areal_rad… |
|
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… |
|
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… |
|
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]*100… |
|
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 |
