 |
tridging_parameterization2.jl - seaice-experiments - sea ice experiments using … |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
LICENSE |
|
--- |
|
tridging_parameterization2.jl (3366B) |
|
--- |
|
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 |
|
9 |
|
10 d = 0.02 # ice particle thickness [m] |
|
11 L = 200*d # ice-floe length |
|
12 E = 3.5e6 # Young's modulus |
|
13 ν = 0.285 # Poisson's ratio [-] |
|
14 ρ_w = 1e3 # Water density [kg/m^3] |
|
15 g = 9.8 # Gravitational acceleration [m/s^2] |
|
16 |
|
17 |
|
18 ########################################################################… |
|
19 #### Compare ridging compressive stress with Granular.jl parameterization |
|
20 |
|
21 h1 = Float64[] # thickness of ice floe 1 [m] |
|
22 h2 = Float64[] # thickness of ice floe 2 [m] |
|
23 comp_vel = Float64[] # compressive velocity [m/s] |
|
24 N_max = Float64[] # compressive stress at yield point [Pa] |
|
25 τ_max = Float64[] # shear stress at yield point [Pa] |
|
26 t_yield = Float64[] # time of yield failure [t] |
|
27 d_yield = Float64[] # compressive distance at yield failure [m] |
|
28 N_late_avg = Float64[] # averaged post-failure compressive stress [Pa] |
|
29 τ_late_avg = Float64[] # averaged post-failure shear stress [Pa] |
|
30 tensile_strength = Float64[] # tensile strength [Pa] |
|
31 |
|
32 nx1 = 100 |
|
33 nx2 = 100 |
|
34 ny1 = 10 |
|
35 ny2 = 15 |
|
36 h1 = ny1*d*sin(π/3) # correct thickness for triangular packing |
|
37 h2 = ny2*d*sin(π/3) # correct thickness for triangular packing |
|
38 cv = 0.2 |
|
39 sim = Granular.createSimulation() |
|
40 fracture_toughness = 1.9e6 |
|
41 r1 = nx1*d |
|
42 r2 = nx2*d |
|
43 coulomb_friction = 0.4 |
|
44 Granular.addGrainCylindrical!(sim, [0., 0.], r1, h1, |
|
45 fracture_toughness=fracture_toughness, |
|
46 youngs_modulus=E, |
|
47 contact_static_friction=coulomb_friction, |
|
48 contact_dynamic_friction=coulomb_friction, |
|
49 lin_vel=[cv, 0.], |
|
50 fixed=true) |
|
51 Granular.addGrainCylindrical!(sim, [r1 + r2 + 12*d, 0.0], r2, h2, |
|
52 fracture_toughness=fracture_toughness, |
|
53 youngs_modulus=E, |
|
54 contact_static_friction=coulomb_friction, |
|
55 contact_dynamic_friction=coulomb_friction, |
|
56 fixed=true) |
|
57 Granular.setTimeStep!(sim) |
|
58 compressive_distance = 3.0 |
|
59 Granular.setTotalTime!(sim, compressive_distance/cv) |
|
60 Granular.removeSimulationFiles(sim) |
|
61 Granular.setOutputFileInterval!(sim, compressive_distance/cv * (1.0/10.0… |
|
62 dist = Float64[] |
|
63 f_n_parameterization = Float64[] |
|
64 N_parameterization = Float64[] |
|
65 #println("Peak stress by fit: $(our_strength_model(min(h1[1], h2[1]), fi… |
|
66 r12 = 2.0*r1*r2/(r1 + r2) |
|
67 Lz_12 = min(h1, h2) |
|
68 while sim.time < sim.time_total |
|
69 for i=1:1 |
|
70 Granular.run!(sim, single_step=true) |
|
71 end |
|
72 append!(dist, cv*sim.time) |
|
73 append!(f_n_parameterization, -sim.grains[1].force[1]) |
|
74 append!(N_parameterization, -sim.grains[1].force[1]/(r12*Lz_12)) |
|
75 end |
|
76 |
|
77 PyPlot.figure(figsize=(4,3)) |
|
78 A_exp = ny1*d * d |
|
79 # Subtract drag against water from data set |
|
80 #PyPlot.plot(data[1,:].*cv, (data[2,:] - mean(data[2,1:100]))./1e3, "C1-… |
|
81 # label="Two-floe experiment") |
|
82 PyPlot.plot(dist, N_parameterization./1e3, "C2-", label="Plan-view param… |
|
83 PyPlot.xlabel("Compressive distance [m]") |
|
84 PyPlot.ylabel("Compressive stress [kPa]") |
|
85 #PyPlot.legend() |
|
86 PyPlot.tight_layout() |
|
87 PyPlot.savefig("ridging_parameterization2.png") |
|
88 PyPlot.savefig("ridging_parameterization2.pdf") |
