ridging_bulk_simulation.jl (9545B)
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 = 1285e3 43 "--r_min" 44 help = "min. grain radius [m]" 45 arg_type = Float64 46 default = 20.0 47 "--r_max" 48 help = "max. grain radius [m]" 49 arg_type = Float64 50 default = 200.0 51 "--thickness" 52 help = "grain thickness [m]" 53 arg_type = Float64 54 default = 1. 55 "--rotating" 56 help = "allow the grains to rotate" 57 arg_type = Bool 58 default = true 59 "--compressive_velocity" 60 help = "compressive velocity [m/s]" 61 arg_type = Float64 62 default = 0.1 63 "--seed" 64 help = "seed for the pseudo RNG" 65 arg_type = Int 66 default = 1 67 "id" 68 help = "simulation id" 69 required = true 70 end 71 return ArgParse.parse_args(s) 72 end 73 74 function report_args(parsed_args) 75 println("Parsed args:") 76 for (arg,val) in parsed_args 77 println(" $arg => $val") 78 end 79 end 80 81 function run_simulation(id::String, 82 E::Float64, 83 nu::Float64, 84 mu_d::Float64, 85 tensile_strength::Float64, 86 tensile_strength_std_dev::Float64, 87 shear_strength::Float64, 88 fracture_toughness::Float64, 89 r_min::Float64, 90 r_max::Float64, 91 thickness::Float64, 92 rotating::Bool, 93 compressive_velocity::Float64, 94 seed::Int) 95 96 ############################################################################ 97 #### SIMULATION PARAMETERS 98 ############################################################################ 99 100 # Common simulation identifier 101 id_prefix = id 102 103 # Grain mechanical properties 104 youngs_modulus = E # Elastic modulus [Pa] 105 poissons_ratio = nu # Shear-stiffness ratio [-] 106 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 107 108 # Simulation duration [s] 109 t_total = 60.0 * 60.0 110 111 # Temporal interval between output files 112 file_dt = t_total/200 113 114 # Misc parameters 115 water_density = 1000. 116 117 # World size [m] 118 L = [1e3, 1e3, 10.0] 119 120 Random.seed!(seed) 121 122 ############################################################################ 123 #### Create ice floes 124 ############################################################################ 125 sim = Granular.createSimulation("$(id_prefix)") 126 Granular.removeSimulationFiles(sim) 127 128 # Create box edges of ice floes with r_min radius, moving inward 129 r_walls = r_min 130 #= for x in LinRange(0.0, L[1] - r_walls, Int(ceil(L[1]/(r_walls*2)))) =# 131 #= Granular.addGrainCylindrical!(sim, [x, r_min], =# 132 #= r_walls, thickness, fixed=true, =# 133 #= lin_vel=[0.0, compressive_velocity/2.], =# 134 #= verbose=false) =# 135 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =# 136 #= r_walls, thickness, fixed=true, =# 137 #= lin_vel=[0.0, -compressive_velocity/2.], =# 138 #= verbose=false) =# 139 #= end =# 140 for y in LinRange(r_walls, L[2] - r_walls, Int(ceil(L[2]/(r_walls*2)))) 141 Granular.addGrainCylindrical!(sim, [r_min, y], 142 r_walls, thickness, fixed=true, 143 #lin_vel=[compressive_velocity/2., 0.0], 144 verbose=false) 145 Granular.addGrainCylindrical!(sim, [L[1] - r_min, y], 146 r_walls, thickness, fixed=true, 147 lin_vel=[-compressive_velocity/2., 0.0], 148 verbose=false) 149 end 150 151 # Create a grid for contact searching spanning the extent of the grains 152 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2] 153 sim.ocean = Granular.createRegularOceanGrid(n, L) 154 Granular.setGridBoundaryConditions!(sim.ocean, "inactive", "-x +x") 155 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "-y +y") 156 157 # Add unconstrained ice floes 158 #= Granular.regularPacking!(sim, =# 159 #= n, =# 160 #= r_min, r_max, =# 161 #= seed=seed) =# 162 Granular.irregularPacking!(sim, 163 radius_max=r_max, radius_min=r_min, 164 thickness=thickness, 165 #binary_radius_search=true, 166 seed=seed) 167 168 # Set grain mechanical properties 169 for grain in sim.grains 170 grain.youngs_modulus = youngs_modulus 171 grain.poissons_ratio = poissons_ratio 172 grain.tensile_strength = abs(tensile_strength + 173 randn()*tensile_strength_std_dev) 174 grain.shear_strength = abs(shear_strength + 175 randn()*tensile_strength_std_dev) 176 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 177 grain.contact_dynamic_friction = contact_dynamic_friction 178 grain.fracture_toughness = fracture_toughness 179 grain.rotating = rotating 180 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 181 end 182 183 # Create GSD plot to signify that simulation is complete 184 Granular.writeSimulation(sim) 185 render && Granular.plotGrainSizeDistribution(sim) 186 187 Granular.setTimeStep!(sim) 188 Granular.setTotalTime!(sim, t_total) 189 Granular.setOutputFileInterval!(sim, file_dt) 190 render && Granular.plotGrains(sim, verbose=true) 191 192 193 ############################################################################ 194 #### RUN THE SIMULATION 195 ############################################################################ 196 #= theta = 0.0; submerged_volume = 0.0 =# 197 time = Float64[] 198 N = Float64[] # normal stress on leftmost fixed particles 199 τ = Float64[] # shear stress on leftmost fixed particles 200 A = sim.grains[1].thickness * L[2] 201 while sim.time < sim.time_total 202 203 append!(time, sim.time) 204 205 # Record cumulative force on fixed particles at the +x boundary 206 normal_force = 0.0 207 tangential_force = 0.0 208 for grain in sim.grains 209 if grain.fixed && grain.lin_pos[1] > 0.2*L[1] 210 normal_force += grain.force[1] 211 tangential_force += grain.force[2] 212 end 213 end 214 append!(N, normal_force/A) # compressive = positive 215 append!(τ, tangential_force/A) 216 217 # Run for 100 time steps before measuring bulk forces 218 for i=1:100 219 Granular.run!(sim, single_step=true) 220 end 221 222 end 223 224 # Plot normal stress and shear stress vs time 225 DelimitedFiles.writedlm(id_prefix * "-data.txt", [time, N, τ]) 226 PyPlot.subplot(211) 227 PyPlot.subplots_adjust(hspace=0.0) 228 ax1 = PyPlot.gca() 229 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels 230 PyPlot.plot(time, N/1e3) 231 PyPlot.ylabel("Compressive stress [kPa]") 232 PyPlot.subplot(212, sharex=ax1) 233 PyPlot.plot(time, τ/1e3) 234 PyPlot.xlabel("Time [s]") 235 PyPlot.ylabel("Shear stress [kPa]") 236 PyPlot.tight_layout() 237 PyPlot.savefig(sim.id * "-time_vs_stress.pdf") 238 PyPlot.clf() 239 240 render && Granular.render(sim, trim=false, animation=true) 241 end 242 243 function main() 244 parsed_args = parse_command_line() 245 report_args(parsed_args) 246 247 seed = parsed_args["seed"] 248 run_simulation(parsed_args["id"] * "-seed" * string(seed), 249 parsed_args["E"], 250 parsed_args["nu"], 251 parsed_args["mu_d"], 252 parsed_args["tensile_strength"], 253 parsed_args["tensile_strength_std_dev"], 254 parsed_args["shear_strength"], 255 parsed_args["fracture_toughness"], 256 parsed_args["r_min"], 257 parsed_args["r_max"], 258 parsed_args["thickness"], 259 parsed_args["rotating"], 260 parsed_args["compressive_velocity"], 261 seed) 262 end 263 264 main() 265 end