simulation_cons.jl (7577B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import JLD 5 import PyPlot 6 import ArgParse 7 8 const render = false 9 10 function parse_command_line() 11 s = ArgParse.ArgParseSettings() 12 ArgParse.@add_arg_table s begin 13 "--E" 14 help = "Young's modulus [Pa]" 15 arg_type = Float64 16 default = 2e7 17 "--nu" 18 help = "Poisson's ratio [-]" 19 arg_type = Float64 20 default = 0.285 21 "--mu_d" 22 help = "dynamic friction coefficient [-]" 23 arg_type = Float64 24 default = 0.3 25 "--gamma_n" 26 help = "contact-normal viscoity [-]" 27 arg_type = Float64 28 default = 0.0 29 "--tensile_strength" 30 help = "the maximum tensile strength [Pa]" 31 arg_type = Float64 32 default = 0. 33 "--r_min" 34 help = "minimum grain radius [m]" 35 arg_type = Float64 36 default = 5. 37 "--r_max" 38 help = "maximum grain radius [m]" 39 arg_type = Float64 40 default = 1.35e3 41 default = 50. 42 "--thickness" 43 help = "grain thickness [m]" 44 arg_type = Float64 45 default = 1. 46 "--rotating" 47 help = "allow the grains to rotate" 48 arg_type = Bool 49 default = true 50 "--shearvel" 51 help = "shear velocity [m/s]" 52 arg_type = Float64 53 default = 1.0 54 "--seed" 55 help = "seed for the pseudo RNG" 56 arg_type = Int 57 default = 1 58 "id" 59 help = "simulation id" 60 required = true 61 end 62 return ArgParse.parse_args(s) 63 end 64 65 function report_args(parsed_args) 66 println("Parsed args:") 67 for (arg,val) in parsed_args 68 println(" $arg => $val") 69 end 70 end 71 72 function run_simulation(id::String, 73 E::Float64, 74 nu::Float64, 75 mu_d::Float64, 76 gamma_n::Float64, 77 tensile_strength::Float64, 78 r_min::Float64, 79 r_max::Float64, 80 thickness::Float64, 81 rotating::Bool, 82 shearvel::Float64, 83 seed::Int) 84 85 ############################################################################ 86 #### SIMULATION PARAMETERS # 87 ############################################################################ 88 89 # Common simulation identifier 90 const id_prefix = id 91 92 # Normal stresses for consolidation and shear [Pa] 93 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3] 94 95 # Simulation duration of individual steps [s] 96 const t_cons = 400.0 97 98 99 ############################################################################ 100 #### Step 2: Consolidate the previous output under a constant normal stress# 101 ############################################################################ 102 tau_p = Float64[] 103 tau_u = Float64[] 104 105 for N in N_list 106 sim_init = Granular.createSimulation("$(id_prefix)-init") 107 sim = Granular.readSimulation(sim_init) 108 109 sim.id = "$(id_prefix)-cons-N$(N)Pa" 110 Granular.removeSimulationFiles(sim) 111 112 # Set all linear and rotational velocities to zero 113 Granular.zeroKinematics!(sim) 114 115 # Add a dynamic wall to the top which adds a normal stress downwards. 116 # The normal of this wall is downwards, and we place it at the top of 117 # the granular assemblage. Here, the fluid drag is also disabled 118 y_top = -Inf 119 for grain in sim.grains 120 grain.contact_viscosity_normal = gamma_n 121 if y_top < grain.lin_pos[2] + grain.contact_radius 122 y_top = grain.lin_pos[2] + grain.contact_radius 123 end 124 end 125 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top, 126 bc="normal stress", 127 normal_stress=-N) 128 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass 129 sim.walls[1].contact_viscosity_normal = 10.0*gamma_n 130 info("Placing top wall at y=$y_top") 131 132 # Resize the grid to span the current state 133 Granular.fitGridToGrains!(sim, sim.ocean) 134 135 # Lock the grains at the very bottom so that the lower boundary is rough 136 y_bot = Inf 137 for grain in sim.grains 138 if y_bot > grain.lin_pos[2] - grain.contact_radius 139 y_bot = grain.lin_pos[2] - grain.contact_radius 140 end 141 end 142 const fixed_thickness = 2.0*r_max 143 for grain in sim.grains 144 if grain.lin_pos[2] <= fixed_thickness 145 grain.fixed = true # set x and y acceleration to zero 146 grain.color = 1 147 end 148 end 149 150 # Set current time to zero and reset output file counter 151 Granular.resetTime!(sim) 152 153 # Set the simulation time to run the consolidation for 154 Granular.setTotalTime!(sim, t_cons) 155 156 # Run the consolidation experiment, and monitor top wall position over 157 # time 158 time = Float64[] 159 compaction = Float64[] 160 effective_normal_stress = Float64[] 161 while sim.time < sim.time_total 162 163 for i=1:100 164 Granular.run!(sim, single_step=true) 165 end 166 167 append!(time, sim.time) 168 append!(compaction, sim.walls[1].pos) 169 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 170 171 end 172 defined_normal_stress = ones(length(effective_normal_stress)) * 173 Granular.getWallNormalStress(sim, stress_type="effective") 174 PyPlot.subplot(211) 175 PyPlot.subplots_adjust(hspace=0.0) 176 ax1 = PyPlot.gca() 177 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 178 PyPlot.plot(time, compaction) 179 PyPlot.ylabel("Top wall height [m]") 180 PyPlot.subplot(212, sharex=ax1) 181 PyPlot.plot(time, defined_normal_stress) 182 PyPlot.plot(time, effective_normal_stress) 183 PyPlot.xlabel("Time [s]") 184 PyPlot.ylabel("Normal stress [Pa]") 185 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 186 PyPlot.clf() 187 writedlm(sim.id * "-data.txt", [time, compaction, 188 effective_normal_stress, 189 defined_normal_stress]) 190 191 # Try to render the simulation if `pvpython` is installed on the system 192 if render 193 Granular.render(sim, trim=false) 194 end 195 196 # Save the simulation state to disk in case we need to reuse the 197 # consolidated state (e.g. different shear velocities below) 198 Granular.writeSimulation(sim) 199 200 end 201 end 202 203 function main() 204 parsed_args = parse_command_line() 205 report_args(parsed_args) 206 207 seed = parsed_args["seed"] 208 run_simulation(parsed_args["id"] * "-seed" * string(seed), 209 parsed_args["E"], 210 parsed_args["nu"], 211 parsed_args["mu_d"], 212 parsed_args["gamma_n"], 213 parsed_args["tensile_strength"], 214 parsed_args["r_min"], 215 parsed_args["r_max"], 216 parsed_args["thickness"], 217 parsed_args["rotating"], 218 parsed_args["shearvel"], 219 seed) 220 end 221 222 main()