compress.jl (8221B)
1 #/usr/bin/env julia 2 ENV["MPLBACKEND"] = "Agg" 3 import Granular 4 import PyPlot 5 import ArgParse 6 using Random 7 import DelimitedFiles 8 9 let 10 render = false 11 12 function parse_command_line() 13 s = ArgParse.ArgParseSettings() 14 ArgParse.@add_arg_table s begin 15 "--E" 16 help = "Young's modulus [Pa]" 17 arg_type = Float64 18 default = 2e7 19 "--nu" 20 help = "Poisson's ratio [-]" 21 arg_type = Float64 22 default = 0.285 23 "--mu_d" 24 help = "dynamic friction coefficient [-]" 25 arg_type = Float64 26 default = 0.3 27 "--tensile_strength" 28 help = "the maximum tensile strength [Pa]" 29 arg_type = Float64 30 default = 400e3 31 "--tensile_strength_std_dev" 32 help = "standard deviation of the maximum tensile strength [Pa]" 33 arg_type = Float64 34 default = 0.0 35 "--shear_strength" 36 help = "the maximum shear strength [Pa]" 37 arg_type = Float64 38 default = 200e3 39 "--fracture_toughness" 40 help = "fracture toughness [Pa*m^{1/2}]" 41 arg_type = Float64 42 default = 1e20 43 #default = 1285e3 44 "--r_min" 45 help = "min. grain radius [m]" 46 arg_type = Float64 47 default = 20.0 48 "--r_max" 49 help = "max. grain radius [m]" 50 arg_type = Float64 51 default = 80.0 52 "--thickness" 53 help = "grain thickness [m]" 54 arg_type = Float64 55 default = 1. 56 "--rotating" 57 help = "allow the grains to rotate" 58 arg_type = Bool 59 default = true 60 "--compressive_velocity" 61 help = "compressive velocity [m/s]" 62 arg_type = Float64 63 default = 0.1 64 "--seed" 65 help = "seed for the pseudo RNG" 66 arg_type = Int 67 default = 1 68 "id" 69 help = "simulation id" 70 required = true 71 end 72 return ArgParse.parse_args(s) 73 end 74 75 function report_args(parsed_args) 76 println("Parsed args:") 77 for (arg,val) in parsed_args 78 println(" $arg => $val") 79 end 80 end 81 82 function run_simulation(id::String, 83 E::Float64, 84 nu::Float64, 85 mu_d::Float64, 86 tensile_strength::Float64, 87 tensile_strength_std_dev::Float64, 88 shear_strength::Float64, 89 fracture_toughness::Float64, 90 r_min::Float64, 91 r_max::Float64, 92 thickness::Float64, 93 rotating::Bool, 94 compressive_velocity::Float64, 95 seed::Int) 96 97 ############################################################################ 98 #### SIMULATION PARAMETERS 99 ############################################################################ 100 101 # Grain mechanical properties 102 youngs_modulus = E # Elastic modulus [Pa] 103 poissons_ratio = nu # Shear-stiffness ratio [-] 104 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 105 106 # Simulation duration [s] 107 t_total = 60.0 * 60.0 * 2.0 108 109 # Temporal interval between output files 110 file_dt = 18.0 111 112 ############################################################################ 113 #### Read prior state 114 ############################################################################ 115 sim = Granular.readSimulation("rotating01-seed1/rotating01-seed1.681.jld2") 116 sim.id = id 117 Granular.resetTime!(sim) 118 Granular.zeroKinematics!(sim) 119 Granular.removeSimulationFiles(sim) 120 121 # Set grain mechanical properties 122 for grain in sim.grains 123 grain.youngs_modulus = youngs_modulus 124 grain.poissons_ratio = poissons_ratio 125 grain.tensile_strength = abs(tensile_strength + 126 randn()*tensile_strength_std_dev) 127 grain.shear_strength = abs(shear_strength + 128 randn()*tensile_strength_std_dev) 129 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 130 grain.contact_dynamic_friction = contact_dynamic_friction 131 grain.fracture_toughness = fracture_toughness 132 grain.rotating = rotating 133 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 134 end 135 136 x_max = -Inf 137 x_min = +Inf 138 y_max = -Inf 139 y_min = +Inf 140 r_max = 0.0 141 for grain in sim.grains 142 r = grain.contact_radius 143 144 if x_min > grain.lin_pos[1] - r 145 x_min = grain.lin_pos[1] - r 146 end 147 148 if x_max < grain.lin_pos[1] + r 149 x_max = grain.lin_pos[1] + r 150 end 151 152 if y_min > grain.lin_pos[2] - r 153 y_min = grain.lin_pos[2] - r 154 end 155 156 if y_max < grain.lin_pos[2] + r 157 y_max = grain.lin_pos[2] + r 158 end 159 160 if r_max < r 161 r_max = r 162 end 163 end 164 165 dx = r_max * 2.0 166 L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx] 167 n = convert(Vector{Int}, floor.(L./dx)) 168 sim.ocean = Granular.createRegularOceanGrid(n, L, 169 origo=[x_min - L[1]/3.0, y_min], 170 time=[0.], name="resized") 171 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x") 172 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y") 173 174 sim.walls[1].bc = "velocity" 175 sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total 176 177 # fix grains adjacent to walls to create no-slip BC 178 fixed_thickness = 2.0 * r_max 179 for grain in sim.grains 180 if grain.lin_pos[2] <= y_min + fixed_thickness 181 grain.fixed = true 182 grain.color = 1 183 elseif grain.lin_pos[2] >= y_max - fixed_thickness 184 grain.allow_y_acc = true 185 grain.fixed = true 186 grain.color = 1 187 end 188 end 189 190 Granular.setTimeStep!(sim) 191 Granular.setTotalTime!(sim, t_total) 192 Granular.setOutputFileInterval!(sim, file_dt) 193 render && Granular.plotGrains(sim, verbose=true) 194 195 ############################################################################ 196 #### RUN THE SIMULATION 197 ############################################################################ 198 199 # Run the consolidation experiment, and monitor top wall position over 200 # time 201 time = Float64[] 202 compaction = Float64[] 203 effective_normal_stress = Float64[] 204 while sim.time < sim.time_total 205 206 for i=1:100 207 Granular.run!(sim, single_step=true) 208 end 209 210 append!(time, sim.time) 211 append!(compaction, sim.walls[1].pos) 212 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 213 214 end 215 Granular.writeSimulation(sim) 216 217 defined_normal_stress = ones(length(effective_normal_stress)) * 218 Granular.getWallNormalStress(sim, stress_type="effective") 219 PyPlot.subplot(211) 220 PyPlot.subplots_adjust(hspace=0.0) 221 ax1 = PyPlot.gca() 222 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 223 PyPlot.plot(time, compaction) 224 PyPlot.ylabel("Top wall height [m]") 225 PyPlot.subplot(212, sharex=ax1) 226 PyPlot.plot(time, defined_normal_stress) 227 PyPlot.plot(time, effective_normal_stress) 228 PyPlot.xlabel("Time [s]") 229 PyPlot.ylabel("Normal stress [Pa]") 230 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 231 PyPlot.clf() 232 DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction, 233 effective_normal_stress, 234 defined_normal_stress]) 235 236 render && Granular.render(sim, trim=false, animation=true) 237 end 238 239 function main() 240 parsed_args = parse_command_line() 241 report_args(parsed_args) 242 243 seed = parsed_args["seed"] 244 run_simulation(parsed_args["id"] * "-seed" * string(seed), 245 parsed_args["E"], 246 parsed_args["nu"], 247 parsed_args["mu_d"], 248 parsed_args["tensile_strength"], 249 parsed_args["tensile_strength_std_dev"], 250 parsed_args["shear_strength"], 251 parsed_args["fracture_toughness"], 252 parsed_args["r_min"], 253 parsed_args["r_max"], 254 parsed_args["thickness"], 255 parsed_args["rotating"], 256 parsed_args["compressive_velocity"], 257 seed) 258 end 259 260 main() 261 end