 |
tridging_plots.jl - seaice-experiments - sea ice experiments using Granular.jl |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
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 wh… |
|
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: ti… |
|
84 |
|
85 n_yw = Int(round(size(data)[2]/4)) # yield detected within first qua… |
|
86 N_max_, t_yield_ = findmax(data[2, 1:n_yw]) |
|
87 |
|
88 # Store scalar values for simulation parameters and resultant stress… |
|
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(… |
|
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, τ_l… |
|
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_amund… |
|
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… |
|
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… |
|
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… |
|
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… |
|
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]))… |
|
589 label="Two-floe experiment") |
|
590 PyPlot.plot(dist, N_parameterization./1e3, "C2--", label="Plan-view para… |
|
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-f… |
|
601 PyPlot.semilogy(dist, abs.(N_parameterization)./1e3, "C2--", label="Plan… |
|
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") |
