tsimulation_mohr-coulomb_plot.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
---
tsimulation_mohr-coulomb_plot.jl (4045B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import PyPlot
4 import LsqFit
5
6 const id_prefix = "mohr_coulomb"
7 const vel_shear = 1.0
8
9 # Normal stresses for consolidation and shear [Pa]
10 #const N_list = [20e3, 40e3, 80e3, 160e3]
11 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3]
12
13 #time = Float64[]
14 #shear_strain = Float64[]
15 shear_stress = Float64[]
16 #dilation = Float64[]
17 tau_p = Float64[]
18 tau_u = Float64[]
19 sim_id = ""
20
21 for kind in ["mu0.3_sigma_c0kPa", "mu0.0_sigma_c200kPa"]
22 for N in N_list
23
24 #mohr_coulomb_mu0.3_sigma_c0kPa.pdf-seed1-shear-N20000.0Pa-vel_shear1.0m-s-data
25 sim_id = "$(id_prefix)_$(kind).pdf-seed1-shear-N$(N)Pa-vel_shear$(vel_shear)m-s"
26 info(sim_id)
27 data = readdlm(sim_id * "-data.txt")
28
29 #time = data[1,:]
30 #shear_strain = data[2,:]
31 shear_stress = data[3,:]
32 #dilation = data[4,:]
33
34 append!(tau_p, maximum(shear_stress[1:100]))
35 append!(tau_u, mean(shear_stress[end-100:end]))
36 end
37 end
38
39 # Plot normal stress vs. peak shear friction
40 PyPlot.figure(figsize=(5, 4))
41 PyPlot.plot(N_list, tau_p[1:length(N_list)]./N_list,
42 "xk",
43 label="Frictional DEM")
44 PyPlot.plot(N_list, tau_p[length(N_list)+1:end]./N_list,
45 "+b",
46 label="Cohesive DEM")
47 PyPlot.ylabel("Peak shear friction \$\\mu_p\$ [-]")
48
49 PyPlot.xlabel("Normal stress \$N\$ [Pa]")
50 PyPlot.xlim(0., maximum(N_list)*1.1)
51 PyPlot.legend()
52 PyPlot.savefig(id_prefix * "-mu_p.pdf")
53 PyPlot.clf()
54
55 # Plot normal stress vs. effective shear friction
56 PyPlot.plot(N_list./1e3, tau_u[1:length(N_list)]./N_list,
57 "xk",
58 label="Frictional DEM")
59 PyPlot.plot(N_list./1e3, tau_u[length(N_list)+1:end]./N_list,
60 "+b",
61 label="Cohesive DEM")
62 PyPlot.xlabel("Normal stress \$N\$ [kPa]")
63 PyPlot.ylabel("Effective shear friction \$\\tau_u/N\$ [-]")
64 #PyPlot.xlim(0., maximum(N_list)*1.1/1e3)
65 #PyPlot.xlim(minimum(N_list)*0.9/1e3, maximum(N_list)*1.1/1e3)
66 PyPlot.legend()
67 PyPlot.title("(b)")
68 PyPlot.savefig(id_prefix * "-mu_eff.pdf")
69 PyPlot.clf()
70
71 ### Plot normal stress vs. ultimate shear stress
72 mohr_coulomb(N, param) = param[2] + param[1]*N
73
74 # Frictional DEM
75 fit = LsqFit.curve_fit(mohr_coulomb, N_list, tau_u[1:length(N_list)], [0.3, 0.0])
76 println("### Frictional DEM ###")
77 println(" μ_u = $(fit.param[1]), C = $(fit.param[2]/1e3) kPa")
78 errors = LsqFit.estimate_errors(fit, 0.95)
79 println(" 95% CI errors:")
80 println(" μ_u: ± $(errors[1])")
81 println(" C: ± $(errors[2]/1e3) kPa")
82 label = @sprintf("\$\\tau_u\$(\$N\$, \$\\mu_u\$ = %.2g±%.2g, \$C\$ = %.2g±%.2g kPa)",
83 fit.param[1], errors[1],
84 fit.param[2]/1e3, errors[2]/1e3)
85 PyPlot.plot(N_list./1e3, tau_u[1:length(N_list)]./1e3,
86 "xk",
87 label="Frictional DEM")
88 PyPlot.plot(linspace(N_list[1], N_list[end])./1e3,
89 mohr_coulomb(linspace(N_list[1], N_list[end]), fit.param)./1e3,
90 "-k", alpha=0.5, label=label)
91
92 # Cohesive DEM
93 fit = LsqFit.curve_fit(mohr_coulomb, N_list, tau_u[length(N_list)+1:end], [0.3, 0.0])
94 println("### Cohesive DEM ###")
95 println(" μ_u = $(fit.param[1]), C = $(fit.param[2]/1e3) kPa")
96 errors = LsqFit.estimate_errors(fit, 0.95)
97 println(" 95% CI errors:")
98 println(" μ_u: ± $(errors[1])")
99 println(" C: ± $(errors[2]/1e3) kPa")
100 label = @sprintf("\$\\tau_u\$(\$N\$, \$\\mu_u\$ = %.2g±%.2g, \$C\$ = %.1f±%.2g kPa)",
101 fit.param[1], errors[1],
102 fit.param[2]/1e3, errors[2]/1e3)
103 PyPlot.plot(N_list./1e3, tau_u[length(N_list)+1:end]./1e3,
104 "+b",
105 label="Cohesive DEM")
106 PyPlot.plot(linspace(N_list[1], N_list[end])./1e3,
107 mohr_coulomb(linspace(N_list[1], N_list[end]), fit.param)./1e3,
108 "--b", alpha=0.5, label=label)
109
110 # Annotation
111 PyPlot.xlabel("Normal stress \$N\$ [kPa]")
112 PyPlot.ylabel("Shear stress \$\\tau_u\$ [kPa]")
113 #PyPlot.xlim(minimum(N_list)*0.9/1e3, maximum(N_list)*1.1/1e3)
114 PyPlot.legend()
115 PyPlot.title("(a)")
116 PyPlot.savefig(id_prefix * "-tau_u.pdf")
117 PyPlot.clf()