ridging_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")