simulation_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()