compress.jl (8238B)
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 = 1.6e-4 134 grain.ocean_drag_coeff_vert = 0.14 135 end 136 137 x_max = -Inf 138 x_min = +Inf 139 y_max = -Inf 140 y_min = +Inf 141 r_max = 0.0 142 for grain in sim.grains 143 r = grain.contact_radius 144 145 if x_min > grain.lin_pos[1] - r 146 x_min = grain.lin_pos[1] - r 147 end 148 149 if x_max < grain.lin_pos[1] + r 150 x_max = grain.lin_pos[1] + r 151 end 152 153 if y_min > grain.lin_pos[2] - r 154 y_min = grain.lin_pos[2] - r 155 end 156 157 if y_max < grain.lin_pos[2] + r 158 y_max = grain.lin_pos[2] + r 159 end 160 161 if r_max < r 162 r_max = r 163 end 164 end 165 166 dx = r_max * 2.0 167 L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx] 168 n = convert(Vector{Int}, floor.(L./dx)) 169 sim.ocean = Granular.createRegularOceanGrid(n, L, 170 origo=[x_min - L[1]/3.0, y_min], 171 time=[0.], name="resized") 172 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x") 173 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y") 174 175 sim.walls[1].bc = "velocity" 176 sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total 177 178 # fix grains adjacent to walls to create no-slip BC 179 fixed_thickness = 2.0 * r_max 180 for grain in sim.grains 181 if grain.lin_pos[2] <= y_min + fixed_thickness 182 grain.fixed = true 183 grain.color = 1 184 elseif grain.lin_pos[2] >= y_max - fixed_thickness 185 grain.allow_y_acc = true 186 grain.fixed = true 187 grain.color = 1 188 end 189 end 190 191 Granular.setTimeStep!(sim) 192 Granular.setTotalTime!(sim, t_total) 193 Granular.setOutputFileInterval!(sim, file_dt) 194 render && Granular.plotGrains(sim, verbose=true) 195 196 ############################################################################ 197 #### RUN THE SIMULATION 198 ############################################################################ 199 200 # Run the consolidation experiment, and monitor top wall position over 201 # time 202 time = Float64[] 203 compaction = Float64[] 204 effective_normal_stress = Float64[] 205 while sim.time < sim.time_total 206 207 for i=1:100 208 Granular.run!(sim, single_step=true) 209 end 210 211 append!(time, sim.time) 212 append!(compaction, sim.walls[1].pos) 213 append!(effective_normal_stress, Granular.getWallNormalStress(sim)) 214 215 end 216 Granular.writeSimulation(sim) 217 218 defined_normal_stress = ones(length(effective_normal_stress)) * 219 Granular.getWallNormalStress(sim, stress_type="effective") 220 PyPlot.subplot(211) 221 PyPlot.subplots_adjust(hspace=0.0) 222 ax1 = PyPlot.gca() 223 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 224 PyPlot.plot(time, compaction) 225 PyPlot.ylabel("Top wall height [m]") 226 PyPlot.subplot(212, sharex=ax1) 227 PyPlot.plot(time, defined_normal_stress) 228 PyPlot.plot(time, effective_normal_stress) 229 PyPlot.xlabel("Time [s]") 230 PyPlot.ylabel("Normal stress [Pa]") 231 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf") 232 PyPlot.clf() 233 DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction, 234 effective_normal_stress, 235 defined_normal_stress]) 236 237 render && Granular.render(sim, trim=false, animation=true) 238 end 239 240 function main() 241 parsed_args = parse_command_line() 242 report_args(parsed_args) 243 244 seed = parsed_args["seed"] 245 run_simulation(parsed_args["id"] * "-seed" * string(seed), 246 parsed_args["E"], 247 parsed_args["nu"], 248 parsed_args["mu_d"], 249 parsed_args["tensile_strength"], 250 parsed_args["tensile_strength_std_dev"], 251 parsed_args["shear_strength"], 252 parsed_args["fracture_toughness"], 253 parsed_args["r_min"], 254 parsed_args["r_max"], 255 parsed_args["thickness"], 256 parsed_args["rotating"], 257 parsed_args["compressive_velocity"], 258 seed) 259 end 260 261 main() 262 end