compaction.jl (2506B)
1 import Granular 2 import JLD2 3 import PyPlot 4 import Dates 5 6 sim = Granular.readSimulation("stacked60k.jld2") 7 SimSettings = JLD2.load("SimSettings.jld2") 8 9 N = 20e3 10 SimSettings["N"] = N 11 12 JLD2.save("SimSettings.jld2", SimSettings) 13 14 t_comp = 10.0 #compaction max duration [s] 15 16 sim.id = "compaction-N$(N)Pa" 17 18 Granular.zeroKinematics!(sim) 19 20 y_top = -Inf 21 for grain in sim.grains 22 grain.contact_viscosity_normal = 0 23 if y_top < grain.lin_pos[2] + grain.contact_radius 24 global y_top = grain.lin_pos[2] + grain.contact_radius 25 end 26 end 27 28 Granular.addWallLinearFrictionless!(sim, [0., 1.],y_top, 29 bc="normal stress", 30 normal_stress=-N, 31 contact_viscosity_normal=1e3) 32 33 Granular.fitGridToGrains!(sim,sim.ocean) 34 35 y_bot = Inf 36 for grain in sim.grains 37 if y_bot > grain.lin_pos[2] - grain.contact_radius 38 global y_bot = grain.lin_pos[2] - grain.contact_radius 39 end 40 end 41 fixed_thickness = 2. * SimSettings["r_max"] 42 for grain in sim.grains 43 if grain.lin_pos[2] <= fixed_thickness 44 grain.fixed = true # set x and y acceleration to zero 45 end 46 end 47 48 Granular.resetTime!(sim) 49 50 Granular.setTotalTime!(sim,t_comp) 51 52 time = Float64[] 53 compaction = Float64[] 54 effective_normal_stress = Float64[] 55 56 filen = 1 57 t_start = Dates.now() 58 59 60 while sim.time < sim.time_total 61 62 for i = 1:100 #run for a while before measuring the state of the top wall 63 Granular.run!(sim, single_step=true) 64 end 65 66 #Granular.writeSimulation(sim,filename = "compaction-N$(N)Pa/comp$(filen).jld2") 67 filen = filen+1 68 69 append!(time, sim.time) 70 append!(compaction, sim.walls[1].pos) 71 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 72 73 t_now = Dates.now() 74 dur = Dates.canonicalize(t_now-t_start) 75 print("Time elapsed: ",dur) 76 end 77 78 defined_normal_stress = ones(length(effective_normal_stress)) * 79 Granular.getWallNormalStress(sim, stress_type="effective") 80 81 PyPlot.subplot(211) 82 PyPlot.subplots_adjust(hspace=0.0) 83 ax1 = PyPlot.gca() 84 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels 85 PyPlot.plot(time, compaction) 86 PyPlot.ylabel("Top wall height [m]") 87 PyPlot.subplot(212, sharex=ax1) 88 PyPlot.plot(time, defined_normal_stress) 89 PyPlot.plot(time, effective_normal_stress) 90 PyPlot.xlabel("Time [s]") 91 PyPlot.ylabel("Normal stress [Pa]") 92 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 93 PyPlot.clf() 94 95 Granular.writeSimulation(sim,filename = "comp.jld2")