init-assemblage.jl (8835B)
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 # Common simulation identifier 102 id_prefix = id 103 104 # Grain mechanical properties 105 youngs_modulus = E # Elastic modulus [Pa] 106 poissons_ratio = nu # Shear-stiffness ratio [-] 107 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-] 108 109 # Simulation duration [s] 110 t_total = 60.0 * 60.0 * 4.0 111 112 # Temporal interval between output files 113 file_dt = 18.0 114 115 # Misc parameters 116 water_density = 1000. 117 118 # World size [m] 119 L = [5e3, 2e4, 10.0] 120 121 Random.seed!(seed) 122 123 ############################################################################ 124 #### Create ice floes 125 ############################################################################ 126 sim = Granular.createSimulation("$(id_prefix)") 127 Granular.removeSimulationFiles(sim) 128 129 # Create box edges of ice floes with r_min radius, moving inward 130 r_walls = r_min 131 #= for x in LinRange(0.0, L[1] - r_walls, Int(ceil(L[1]/(r_walls*2)))) =# 132 #= Granular.addGrainCylindrical!(sim, [x, r_min], =# 133 #= r_walls, thickness, fixed=true, =# 134 #= lin_vel=[0.0, compressive_velocity/2.], =# 135 #= verbose=false) =# 136 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =# 137 #= r_walls, thickness, fixed=true, =# 138 #= lin_vel=[0.0, -compressive_velocity/2.], =# 139 #= verbose=false) =# 140 #= end =# 141 #for y in LinRange(r_walls, L[2] - r_walls, Int(ceil(L[2]/(r_walls*2)))) 142 # Granular.addGrainCylindrical!(sim, [r_min, y], 143 # r_walls, thickness, fixed=true, 144 # #lin_vel=[compressive_velocity/2., 0.0], 145 # verbose=false) 146 # Granular.addGrainCylindrical!(sim, [L[1] - r_min, y], 147 # r_walls, thickness, fixed=true, 148 # lin_vel=[-compressive_velocity/2., 0.0], 149 # verbose=false) 150 #end 151 152 # Create a grid for contact searching spanning the extent of the grains 153 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2] 154 sim.ocean = Granular.createRegularOceanGrid(n, L) 155 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x") 156 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y") 157 158 # Add unconstrained ice floes 159 #= Granular.regularPacking!(sim, =# 160 #= n, =# 161 #= r_min, r_max, =# 162 #= seed=seed) =# 163 Granular.irregularPacking!(sim, 164 radius_max=r_max, radius_min=r_min, 165 thickness=thickness, 166 #binary_radius_search=true, 167 seed=seed) 168 169 # Set grain mechanical properties 170 for grain in sim.grains 171 grain.youngs_modulus = youngs_modulus 172 grain.poissons_ratio = poissons_ratio 173 grain.tensile_strength = abs(tensile_strength + 174 randn()*tensile_strength_std_dev) 175 grain.shear_strength = abs(shear_strength + 176 randn()*tensile_strength_std_dev) 177 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 178 grain.contact_dynamic_friction = contact_dynamic_friction 179 grain.fracture_toughness = fracture_toughness 180 grain.rotating = rotating 181 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 182 end 183 184 # Create GSD plot to signify that simulation is complete 185 Granular.writeSimulation(sim) 186 render && Granular.plotGrainSizeDistribution(sim) 187 188 # Add a dynamic wall to the top which adds a normal stress downwards. 189 # The normal of this wall is downwards, and we place it at the top of 190 # the granular assemblage. Here, the fluid drag is also disabled 191 y_top = -Inf 192 for grain in sim.grains 193 if y_top < grain.lin_pos[2] + grain.contact_radius 194 y_top = grain.lin_pos[2] + grain.contact_radius 195 end 196 end 197 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top, 198 bc="normal stress", 199 normal_stress=-20e3) 200 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass 201 @info("Placing top wall at y=$y_top") 202 203 # Resize the grid to span the current state 204 Granular.fitGridToGrains!(sim, sim.ocean) 205 206 Granular.setTimeStep!(sim) 207 Granular.setTotalTime!(sim, t_total) 208 Granular.setOutputFileInterval!(sim, file_dt) 209 render && Granular.plotGrains(sim, verbose=true) 210 211 ############################################################################ 212 #### RUN THE SIMULATION 213 ############################################################################ 214 215 Granular.run!(sim) 216 Granular.writeSimulation(sim) 217 218 render && Granular.render(sim, trim=false, animation=true) 219 end 220 221 function main() 222 parsed_args = parse_command_line() 223 report_args(parsed_args) 224 225 seed = parsed_args["seed"] 226 run_simulation(parsed_args["id"] * "-seed" * string(seed), 227 parsed_args["E"], 228 parsed_args["nu"], 229 parsed_args["mu_d"], 230 parsed_args["tensile_strength"], 231 parsed_args["tensile_strength_std_dev"], 232 parsed_args["shear_strength"], 233 parsed_args["fracture_toughness"], 234 parsed_args["r_min"], 235 parsed_args["r_max"], 236 parsed_args["thickness"], 237 parsed_args["rotating"], 238 parsed_args["compressive_velocity"], 239 seed) 240 end 241 242 main() 243 end