tridging_plots.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_plots.jl (22903B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import PyPlot
4 import LsqFit
5 import Granular
6 import DelimitedFiles
7 import Statistics
8 import LinearAlgebra
9
10 id_prefix = "ridging-seed1"
11 #compressive_velocity = [0.2, 0.1, 0.05]
12 compressive_velocity = [0.1]
13 ny_fixed = 10
14 #ny_variable = [3, 5, 7, 9, 11, 13, 15, 17, 19, 21]
15 ny_variable = [3, 5, 7, 9, 11]
16 # I am missing ridging-seed1-ny1_10-ny2_*-cv_0.1-seed1-data.txt files where ny2 = [3,5,7,9]
17
18 d = 0.02 # ice particle thickness [m]
19 L = 200*d # ice-floe length
20 E = 2e7 # Young's modulus
21 ν = 0.285 # Poisson's ratio [-]
22 ρ_w = 1e3 # Water density [kg/m^3]
23 g = 9.8 # Gravitational acceleration [m/s^2]
24
25 lbl_amundrud = "Elastic sheet,"
26 lbl_amundrud_m1 = "$lbl_amundrud \$m=1\$"
27 lbl_amundrud_m2 = "$lbl_amundrud \$m=2\$"
28 lbl_amundrud_m3 = "$lbl_amundrud \$m=3\$"
29 lbl_amundrud_m4 = "$lbl_amundrud \$m=4\$"
30 lbl_amundrud_m5 = "$lbl_amundrud \$m=5\$"
31
32 lbl_amundrud_adj = "Elastic beam,"
33 lbl_amundrud_adj_m1 = "$lbl_amundrud_adj \$m=1\$"
34 lbl_amundrud_adj_m2 = "$lbl_amundrud_adj \$m=2\$"
35 lbl_amundrud_adj_m3 = "$lbl_amundrud_adj \$m=3\$"
36 lbl_amundrud_adj_m4 = "$lbl_amundrud_adj \$m=4\$"
37 lbl_amundrud_adj_m5 = "$lbl_amundrud_adj \$m=5\$"
38
39 function buckling_strength(E::Float64, ν::Float64, h::Any, L::Float64,
40 ρ_w::Float64, g::Float64; m::Int=1,
41 method::String="Amundrud")
42
43 draft = h .* 0.9
44 if method == "Amundrud" || method == "Hopkins-Adjusted"
45 strength = E .* π^2 * (m^2 + 1)^2 ./ (12 * (1 - ν^2) * m^2) .*
46 draft.^2.0/L^2 + ρ_w * g ./ m^2 .* L^2 ./ draft
47
48 if method == "Hopkins-Adjusted"
49 return strength ./ (4*π^2/(1 - ν^2))^0.25
50 else
51 return strength
52 end
53
54 elseif method == "Amundrud2"
55 return 2 * π * (m^2 + 1 / m^2) .*
56 sqrt.(E * ρ_w * g .* draft ./ (12*(1 - ν^2)))
57
58 elseif method == "Parmerter"
59 return (E * ρ_w * g/(12*(1 - ν^2)))^0.5 .* h.^1.5 # N/m
60
61 elseif method == "Hopkins"
62 return ρ_w * g * sqrt.(E .* h.^3 / (12 * ρ_w * 9.8)) # N/m
63
64 else
65 error("Method $(method) not understood")
66 end
67 end
68
69 function readTimeSeries(id::String,
70 ny1::Int,
71 ny2::Int,
72 cv::Float64,
73 h1::Vector{Float64},
74 h2::Vector{Float64},
75 comp_vel::Vector{Float64},
76 N_max::Vector{Float64},
77 τ_max::Vector{Float64},
78 t_yield::Vector{Float64},
79 d_yield::Vector{Float64},
80 N_late_avg::Vector{Float64},
81 τ_late_avg::Vector{Float64})
82
83 data = DelimitedFiles.readdlm(id * "-seed1-data.txt") # file is: time, N, τ
84
85 n_yw = Int(round(size(data)[2]/4)) # yield detected within first quarter
86 N_max_, t_yield_ = findmax(data[2, 1:n_yw])
87
88 # Store scalar values for simulation parameters and resultant stresses
89 h1_ = ny1*d*sin(π/3) # correct thickness for triangular packing
90 h2_ = ny2*d*sin(π/3) # correct thickness for triangular packing
91 append!(h1, h1_)
92 append!(h2, h2_)
93 append!(comp_vel, cv)
94 append!(N_max, N_max_)
95 append!(τ_max, maximum(abs.(data[3, 1:n_yw])))
96 append!(t_yield, data[1,:][t_yield_])
97 append!(d_yield, t_yield_*cv)
98 append!(N_late_avg, Statistics.mean(data[2, n_yw:end]))
99 append!(τ_late_avg, Statistics.mean(data[3, n_yw:end]))
100
101 return data
102
103 # Plot
104 #PyPlot.plot(data[1,:]*cv, data[2,:]/1e3 / min(h1_, h2_),
105 #PyPlot.plot(data[1,:]*cv, data[2,:]/1e3,
106 PyPlot.plot(data[1,:]*cv, data[2,:]/1e3 .* min(h1_, h2_)^1.5 ./ min(h1_, h2_),
107 linewidth=0.2, alpha=0.5)
108 #color="gray", linewidth=0.2, alpha=0.5)
109 end
110
111 h1 = Float64[] # thickness of ice floe 1 [m]
112 h2 = Float64[] # thickness of ice floe 2 [m]
113 comp_vel = Float64[] # compressive velocity [m/s]
114 N_max = Float64[] # compressive stress at yield point [Pa]
115 τ_max = Float64[] # shear stress at yield point [Pa]
116 t_yield = Float64[] # time of yield failure [t]
117 d_yield = Float64[] # compressive distance at yield failure [m]
118 N_late_avg = Float64[] # averaged post-failure compressive stress [Pa]
119 τ_late_avg = Float64[] # averaged post-failure shear stress [Pa]
120 tensile_strength = Float64[] # tensile strength [Pa]
121
122 PyPlot.figure(figsize=(4,3))
123 for cv in compressive_velocity
124 for ny2 in ny_variable
125 readTimeSeries(id_prefix * "-ny1_$(ny_fixed)-ny2_$(ny2)-cv_$(cv)",
126 ny_fixed,
127 ny2,
128 cv,
129 h1,
130 h2,
131 comp_vel,
132 N_max,
133 τ_max,
134 t_yield,
135 d_yield,
136 N_late_avg,
137 τ_late_avg)
138 append!(tensile_strength, 400e3)
139 end
140 #=for ny1 in ny_variable
141 (ny1 == 21 && cv ≈ 0.1) && continue
142 readTimeSeries(id_prefix * "-ny1_$(ny1)-ny2_$(ny_fixed)-cv_$(cv)",
143 ny1,
144 ny_fixed,
145 cv,
146 h1,
147 h2,
148 comp_vel,
149 N_max,
150 τ_max,
151 t_yield,
152 d_yield,
153 N_late_avg,
154 τ_late_avg)
155 append!(tensile_strength, 400e3)
156 end=#
157 end
158
159 ## Resolution tests: ridging_res-seed1-size_scaling_{0.5,1.0,2.0}-seed1
160 #for size_scaling in [0.5, 1.0, 2.0]
161 for size_scaling in [0.5, 1.0]
162 cv = 0.1
163 d = 0.02/size_scaling
164 ny1 = Int(round(10*size_scaling))
165 ny2 = Int(round(11*size_scaling))
166 readTimeSeries("ridging_res-seed1-size_scaling_$(size_scaling)",
167 ny1,
168 ny2,
169 cv,
170 h1,
171 h2,
172 comp_vel,
173 N_max,
174 τ_max,
175 t_yield,
176 d_yield,
177 N_late_avg,
178 τ_late_avg)
179 append!(tensile_strength, 400e3)
180 end
181
182 PyPlot.legend()
183 PyPlot.xlabel("Compressive distance \$d\$ [m]")
184 PyPlot.ylabel("Compressive stress \$N\$ [kPa]")
185 PyPlot.tight_layout()
186 PyPlot.savefig(id_prefix * "-N.pdf")
187 PyPlot.clf()
188
189 for ny2 in [3, 5, 7, 9, 11]
190 readTimeSeries("elastoridge-seed1-ny1_$(ny_fixed)-ny2_$(ny2)-cv_0.1",
191 ny_fixed,
192 ny2,
193 0.1,
194 h1,
195 h2,
196 comp_vel,
197 N_max,
198 τ_max,
199 t_yield,
200 d_yield,
201 N_late_avg,
202 τ_late_avg)
203 append!(tensile_strength, 1e24)
204 end
205 for ny in [3, 5, 7, 9, 10]
206 readTimeSeries("elastobuckle-seed1-ny1_$(ny)-cv_0.1",
207 ny,
208 ny,
209 0.1,
210 h1,
211 h2,
212 comp_vel,
213 N_max,
214 τ_max,
215 t_yield,
216 d_yield,
217 N_late_avg,
218 τ_late_avg)
219 append!(tensile_strength, 2e24)
220 end
221
222 # Write summary results to text file
223 DelimitedFiles.writedlm(id_prefix * "-data.txt",
224 [h1, h2, comp_vel, N_max, τ_max, t_yield, d_yield, N_late_avg, τ_late_avg])
225 #PyPlot.subplot(211)
226 #PyPlot.subplots_adjust(hspace=0.0)
227 #ax1 = PyPlot.gca()
228 #PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
229
230 PyPlot.figure(figsize=(5,4))
231
232
233 ## N_max
234 for cv in compressive_velocity
235 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
236 PyPlot.plot(h1[I]+h2[I], N_max[I]/1e3, "+",
237 label="\$c_\\mathrm{v}\$ = $cv m/s")
238 end
239 PyPlot.legend()
240 PyPlot.xlabel("Cumulative ice thickness \$h_1 + h_2\$ [m]")
241 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
242 PyPlot.tight_layout()
243 PyPlot.savefig(id_prefix * "-h1+h2-N_max.pdf")
244 PyPlot.clf()
245
246
247
248
249 PyPlot.figure(figsize=(4,4))
250 #fig, ax = PyPlot.subplots(figsize=(4,4))
251 #PyPlot.subplot(212)
252 h_range = LinRange(0.041, 0.18, 50)
253 #N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=1,
254 #method="Amundrud")
255 #PyPlot.plot(h_range, N_amundrud/1e3, "-k", label=lbl_amundrud_m1)
256 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=2,
257 method="Amundrud")
258 PyPlot.plot(h_range, N_amundrud/1e3, "-k", label=lbl_amundrud_m2)
259 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=3,
260 method="Amundrud")
261 PyPlot.plot(h_range, N_amundrud/1e3, "--k", label=lbl_amundrud_m3)
262 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=4,
263 method="Amundrud")
264 PyPlot.plot(h_range, N_amundrud/1e3, "-.k", label=lbl_amundrud_m4)
265 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=5,
266 method="Amundrud")
267 PyPlot.plot(h_range, N_amundrud/1e3, ":k", label=lbl_amundrud_m5)
268
269 #PyPlot.plot(h_range, N_amundrud/1e3, ":", color="gray", label=lbl_amundrud_m5)
270
271 #= N_parmerter = buckling_strength(E, ν, h_range, L, ρ_w, g, m=0, =#
272 #= method="Parmerter") =#
273 #= PyPlot.plot(h_range, N_parmerter/1e3, "-b", label="Parmerter 1974") =#
274
275 #= N_parmerter = buckling_strength(E, ν, h_range, L, ρ_w, g, m=0, =#
276 #= method="Hopkins") =#
277 #= PyPlot.plot(h_range, N_parmerter/1e3, "--y", label="Hopkins 1998") =#
278
279 #= N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=2, =#
280 #= method="Hopkins-Adjusted") =#
281 #= PyPlot.plot(h_range, N_amundrud/1e3, "-g", label=lbl_amundrud_adj_m2) =#
282 #= N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=3, =#
283 #= method="Hopkins-Adjusted") =#
284 #= PyPlot.plot(h_range, N_amundrud/1e3, "--g", label=lbl_amundrud_adj_m3) =#
285 #= N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=4, =#
286 #= method="Hopkins-Adjusted") =#
287 #= PyPlot.plot(h_range, N_amundrud/1e3, "-.g", label=lbl_amundrud_adj_m4) =#
288 #= N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=5, =#
289 #= method="Hopkins-Adjusted") =#
290 #= PyPlot.plot(h_range, N_amundrud/1e3, ":g", label=lbl_amundrud_adj_m5) =#
291
292 I = findall(tensile_strength .≈ 2e24)
293 PyPlot.plot(min.(h1[I],h2[I]), N_max[I]/1e3, "x",
294 label="Elastic sheet")
295 I = findall(tensile_strength .≈ 1e24)
296 PyPlot.plot(min.(h1[I],h2[I]), N_max[I]/1e3, "o",
297 label="Elastic floes", zorder=10)
298 for cv in compressive_velocity
299 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
300 PyPlot.plot(min.(h1[I],h2[I]), N_max[I]/1e3, "+k",
301 label="Elastic-plastic floes",
302 zorder=11)
303 end
304
305 #our_strength_model(h, fracture_toughness) = fracture_toughness[1].*sqrt.(h)
306 # returns max force [N]
307 our_strength_model(h, fracture_toughness) = fracture_toughness[1].*h.^1.5
308 fracture_toughness = [1e4] # initial guess
309
310 fit = LsqFit.curve_fit(our_strength_model,
311 min.(h1[I],h2[I]),
312 N_max[I],
313 #N_max[I].*min.(h1[I],h2[I]).*d, # stress to force
314 fracture_toughness)
315 covar = LsqFit.estimate_covar(fit)
316 std_error = sqrt.(abs.(LinearAlgebra.diag(covar)))
317 sample_std_dev = std_error .* sqrt(length(I))
318 errors = LsqFit.estimate_errors(fit, 0.95)
319
320 if fit.converged
321 # 95% CI range for normal distribution
322 lower_CI = our_strength_model(h_range, fit.param .- 1.96*std_error[1])
323 upper_CI = our_strength_model(h_range, fit.param .+ 1.96*std_error[1])
324 PyPlot.fill_between(h_range, lower_CI./1e3, upper_CI./1e3, facecolor="y",
325 interpolate=true, alpha=0.33, zorder=1)
326
327 println("fitted value: $(fit.param[1])")
328 PyPlot.plot(h_range, our_strength_model(h_range, fit.param)./1e3, "y:",
329 label="\$K_{Ic}\$min(\$h_1,h_2\$)\$^{3/2}A_{1,2}^{-1}\$",
330 zorder=2)
331 end
332
333 #PyPlot.ylim([0, 1400])
334 PyPlot.xlim([h_range[1], h_range[end]])
335 PyPlot.ylim([0, 800])
336 PyPlot.legend()
337 #ax[:legend](fontsize=10, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
338 PyPlot.legend(fontsize=8, bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
339 ncol=2, mode="expand", borderaxespad=0.)
340 #PyPlot.legend()
341 PyPlot.xlabel("Minimum ice thickness, \$\\min(h_1, h_2)\$ [m]")
342 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
343 #PyPlot.tight_layout()
344 PyPlot.tight_layout(rect=[0, 0, 1, 0.75])
345 PyPlot.savefig(id_prefix * "-min_h1_h2-N_max.pdf", bbox_inches="tight")
346 PyPlot.clf()
347
348
349
350
351 for cv in compressive_velocity
352 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
353 PyPlot.plot(max.(h1[I],h2[I]), N_max[I]/1e3, "+",
354 label="Two elastic-plastic, \$c_\\mathrm{v}\$ = $cv m/s")
355 end
356 I = findall(tensile_strength .≈ 1e24)
357 PyPlot.plot(max.(h1[I],h2[I]), N_max[I]/1e3, "o",
358 label="Two elastic floes, \$c_\\mathrm{v}\$ = 0.1 m/s")
359 I = findall(tensile_strength .≈ 2e24)
360 PyPlot.plot(max.(h1[I],h2[I]), N_max[I]/1e3, "x",
361 label="Single elastic sheet, \$c_\\mathrm{v}\$ = 0.1 m/s")
362 h_range = LinRange(0.17, 0.375, 50)
363 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=1,
364 method="Amundrud")
365 PyPlot.plot(h_range, N_amundrud/1e3, "-k", label=lbl_amundrud_m1)
366 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=2,
367 method="Amundrud")
368 PyPlot.plot(h_range, N_amundrud/1e3, "--k", label=lbl_amundrud_m2)
369 N_amundrud = buckling_strength(E, ν, h_range, L, ρ_w, g, m=3,
370 method="Amundrud")
371 PyPlot.plot(h_range, N_amundrud/1e3, ":k", label=lbl_amundrud_m3)
372 PyPlot.legend()
373 PyPlot.xlabel("Maximum ice thickness \$\\max(h_1, h_2)\$ [m]")
374 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
375 PyPlot.tight_layout()
376 PyPlot.savefig(id_prefix * "-max_h1_h2-N_max.pdf")
377 PyPlot.clf()
378
379 for cv in compressive_velocity
380 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
381 PyPlot.plot(h1[I]./h2[I], N_max[I]/1e3, "+",
382 label="\$c_\\mathrm{v}\$ = $cv m/s")
383 end
384 PyPlot.legend()
385 PyPlot.xlabel("Ice thickness ratio \$h_1/h_2\$ [-]")
386 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
387 PyPlot.tight_layout()
388 PyPlot.savefig(id_prefix * "-max_h1_div_h2-N_max.pdf")
389 PyPlot.clf()
390
391 for cv in compressive_velocity
392 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
393 PyPlot.plot(h1[I].*h2[I], N_max[I]/1e3, "+",
394 label="\$c_\\mathrm{v}\$ = $cv m/s")
395 end
396 PyPlot.legend()
397 PyPlot.xlabel("Ice thickness product \$h_1 \\times h_2\$ [m\$^2\$]")
398 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
399 PyPlot.tight_layout()
400 PyPlot.savefig(id_prefix * "-max_h1_times_h2-N_max.pdf")
401 PyPlot.clf()
402
403 ## N_late_avg
404 for cv in compressive_velocity
405 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
406 PyPlot.plot(h1[I]+h2[I], N_late_avg[I]/1e3, "+",
407 label="\$c_\\mathrm{v}\$ = $cv m/s")
408 end
409 PyPlot.legend()
410 PyPlot.xlabel("Cumulative ice thickness \$h_1 + h_2\$ [m]")
411 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
412 PyPlot.tight_layout()
413 PyPlot.savefig(id_prefix * "-h1+h2-N_late_avg.pdf")
414 PyPlot.clf()
415
416 for cv in compressive_velocity
417 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
418 PyPlot.plot(min.(h1[I],h2[I]), N_late_avg[I]/1e3, "+",
419 label="\$c_\\mathrm{v}\$ = $cv m/s")
420 end
421 PyPlot.legend()
422 PyPlot.xlabel("Minimum ice thickness \$\\min(h_1, h_2)\$ [m]")
423 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
424 PyPlot.tight_layout()
425 PyPlot.savefig(id_prefix * "-min_h1_h2-N_late_avg.pdf")
426 PyPlot.clf()
427
428 for cv in compressive_velocity
429 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
430 PyPlot.plot(max.(h1[I],h2[I]), N_late_avg[I]/1e3, "+",
431 label="\$c_\\mathrm{v}\$ = $cv m/s")
432 end
433 PyPlot.legend()
434 PyPlot.xlabel("Maximum ice thickness \$\\max(h_1, h_2)\$ [m]")
435 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
436 PyPlot.tight_layout()
437 PyPlot.savefig(id_prefix * "-max_h1_h2-N_late_avg.pdf")
438 PyPlot.clf()
439
440 for cv in compressive_velocity
441 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
442 PyPlot.plot(h1[I]./h2[I], N_late_avg[I]/1e3, "+",
443 label="\$c_\\mathrm{v}\$ = $cv m/s")
444 end
445 PyPlot.legend()
446 PyPlot.xlabel("Ice thickness ratio \$h_1/h_2\$ [-]")
447 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
448 PyPlot.tight_layout()
449 PyPlot.savefig(id_prefix * "-max_h1_div_h2-N_late_avg.pdf")
450 PyPlot.clf()
451
452 for cv in compressive_velocity
453 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
454 PyPlot.plot(h1[I].*h2[I], N_late_avg[I]/1e3, "+",
455 label="\$c_\\mathrm{v}\$ = $cv m/s")
456 end
457 PyPlot.legend()
458 PyPlot.xlabel("Ice thickness product \$h_1 \\times h_2\$ [m\$^2\$]")
459 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
460 PyPlot.tight_layout()
461 PyPlot.savefig(id_prefix * "-max_h1_times_h2-N_late_avg.pdf")
462 PyPlot.clf()
463
464 ## comp_vel
465 for cv in compressive_velocity
466 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
467 PyPlot.plot(comp_vel[I], N_max[I]/1e3, "+",
468 label="\$c_\\mathrm{v}\$ = $cv m/s")
469 end
470 PyPlot.legend()
471 PyPlot.xlabel("Compressive velocity \$v_\\mathrm{c}\$ [m/s]")
472 PyPlot.ylabel("Avg. post-failure compressive stress \$\\bar{N}\$ [kPa]")
473 PyPlot.tight_layout()
474 PyPlot.savefig(id_prefix * "-cv-N_max.pdf")
475 PyPlot.clf()
476
477 for cv in compressive_velocity
478 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
479 PyPlot.plot(comp_vel[I], N_late_avg[I]/1e3, "+",
480 label="\$c_\\mathrm{v}\$ = $cv m/s")
481 end
482 PyPlot.legend()
483 PyPlot.xlabel("Compressive velocity \$v_\\mathrm{c}\$ [m/s]")
484 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
485 PyPlot.tight_layout()
486 PyPlot.savefig(id_prefix * "-cv-N_late_avg.pdf")
487 PyPlot.clf()
488
489 ## comp_vel
490 for cv in compressive_velocity
491 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
492 PyPlot.plot(comp_vel[I], N_max[I]/1e3, "+",
493 label="\$c_\\mathrm{v}\$ = $cv m/s")
494 end
495 PyPlot.legend()
496 PyPlot.xlabel("Compressive velocity \$v_\\mathrm{c}\$ [m/s]")
497 PyPlot.ylabel("Max. compressive stress \$N_\\mathrm{max}\$ [kPa]")
498 PyPlot.tight_layout()
499 PyPlot.savefig(id_prefix * "-cv-N_max.pdf")
500 PyPlot.clf()
501
502 for cv in compressive_velocity
503 I = findall((comp_vel .≈ cv) .& (tensile_strength .≈ 400e3))
504 PyPlot.plot(comp_vel[I], N_late_avg[I]/1e3, "+",
505 label="\$c_\\mathrm{v}\$ = $cv m/s")
506 end
507 PyPlot.legend()
508 PyPlot.xlabel("Compressive velocity \$v_\\mathrm{c}\$ [m/s]")
509 PyPlot.ylabel("Avg. post-failure shear stress \$\\bar{N}\$ [kPa]")
510 PyPlot.tight_layout()
511 PyPlot.savefig(id_prefix * "-cv-N_late_avg.pdf")
512 PyPlot.clf()
513
514
515 ################################################################################
516 #### Compare ridging compressive stress with Granular.jl parameterization
517
518 h1 = Float64[] # thickness of ice floe 1 [m]
519 h2 = Float64[] # thickness of ice floe 2 [m]
520 comp_vel = Float64[] # compressive velocity [m/s]
521 N_max = Float64[] # compressive stress at yield point [Pa]
522 τ_max = Float64[] # shear stress at yield point [Pa]
523 t_yield = Float64[] # time of yield failure [t]
524 d_yield = Float64[] # compressive distance at yield failure [m]
525 N_late_avg = Float64[] # averaged post-failure compressive stress [Pa]
526 τ_late_avg = Float64[] # averaged post-failure shear stress [Pa]
527 tensile_strength = Float64[] # tensile strength [Pa]
528
529 nx1 = 100
530 nx2 = 100
531 ny1 = 7
532 ny2 = 10
533 cv = 0.1
534 data = readTimeSeries(id_prefix * "-ny1_$(ny1)-ny2_$(ny2)-cv_$(cv)",
535 ny1,
536 ny2,
537 cv,
538 h1,
539 h2,
540 comp_vel,
541 N_max,
542 τ_max,
543 t_yield,
544 d_yield,
545 N_late_avg,
546 τ_late_avg)
547 sim = Granular.createSimulation()
548 fracture_toughness = 1.42e6
549 r1 = nx1*d
550 r2 = nx2*d
551 coulomb_friction = 0.4
552 Granular.addGrainCylindrical!(sim, [0., 0.], r1, h1[1],
553 fracture_toughness=fracture_toughness,
554 youngs_modulus=2e7,
555 contact_static_friction=coulomb_friction,
556 contact_dynamic_friction=coulomb_friction,
557 lin_vel=[cv, 0.],
558 fixed=true)
559 Granular.addGrainCylindrical!(sim, [r1 + r2 + 5*d, 0.0], r2, h2[1],
560 fracture_toughness=fracture_toughness,
561 youngs_modulus=2e7,
562 contact_static_friction=coulomb_friction,
563 contact_dynamic_friction=coulomb_friction,
564 fixed=true)
565 Granular.setTimeStep!(sim)
566 compressive_distance = 2.0
567 Granular.setTotalTime!(sim, compressive_distance/cv)
568 Granular.removeSimulationFiles(sim)
569 Granular.setOutputFileInterval!(sim, compressive_distance/cv * (1/10))
570 dist = Float64[]
571 f_n_parameterization = Float64[]
572 N_parameterization = Float64[]
573 println("Peak stress by fit: $(our_strength_model(min(h1[1], h2[1]), fit.param[1])) Pa")
574 r12 = 2.0*r1*r2/(r1 + r2)
575 Lz_12 = min(h1[1], h2[1])
576 while sim.time < sim.time_total
577 for i=1:1
578 Granular.run!(sim, single_step=true)
579 end
580 append!(dist, cv*sim.time)
581 append!(f_n_parameterization, -sim.grains[1].force[1])
582 append!(N_parameterization, -sim.grains[1].force[1]/(r12*Lz_12))
583 end
584
585 PyPlot.figure(figsize=(4,3))
586 A_exp = ny1*d * d
587 # Subtract drag against water from data set
588 PyPlot.plot(data[1,:].*cv, (data[2,:] .- Statistics.mean(data[2,1:100]))./1e3, "C1-",
589 label="Two-floe experiment")
590 PyPlot.plot(dist, N_parameterization./1e3, "C2--", label="Plan-view parameterization")
591 PyPlot.xlabel("Compressive distance [m]")
592 PyPlot.ylabel("Compressive stress [kPa]")
593 PyPlot.legend()
594 PyPlot.tight_layout()
595 PyPlot.savefig("ridging_parameterization.png")
596 PyPlot.savefig("ridging_parameterization.pdf")
597
598 PyPlot.clf()
599
600 PyPlot.semilogy(data[1,:].*cv, abs.(data[2,:])./1e3, "C1-", label="Two-floe experiment")
601 PyPlot.semilogy(dist, abs.(N_parameterization)./1e3, "C2--", label="Plan-view parameterization")
602 PyPlot.xlabel("Compressive distance [m]")
603 PyPlot.ylabel("Compressive stress [kPa]")
604 PyPlot.legend()
605 PyPlot.tight_layout()
606 PyPlot.savefig("ridging_parameterization_semilog.png")
607 PyPlot.savefig("ridging_parameterization_semilog.pdf")