init-assemblage.jl (7396B)
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 132 # Create a grid for contact searching spanning the extent of the grains 133 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2] 134 sim.ocean = Granular.createRegularOceanGrid(n, L) 135 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x") 136 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y") 137 138 # Add unconstrained ice floes 139 Granular.irregularPacking!(sim, 140 radius_max=r_max, radius_min=r_min, 141 thickness=thickness, 142 #binary_radius_search=true, 143 seed=seed) 144 145 # Set grain mechanical properties 146 for grain in sim.grains 147 grain.youngs_modulus = youngs_modulus 148 grain.poissons_ratio = poissons_ratio 149 grain.tensile_strength = abs(tensile_strength + 150 randn()*tensile_strength_std_dev) 151 grain.shear_strength = abs(shear_strength + 152 randn()*tensile_strength_std_dev) 153 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d 154 grain.contact_dynamic_friction = contact_dynamic_friction 155 grain.fracture_toughness = fracture_toughness 156 grain.rotating = rotating 157 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends 158 end 159 160 # Create GSD plot to signify that simulation is complete 161 Granular.writeSimulation(sim) 162 render && Granular.plotGrainSizeDistribution(sim) 163 164 # Add a dynamic wall to the top which adds a normal stress downwards. 165 # The normal of this wall is downwards, and we place it at the top of 166 # the granular assemblage. Here, the fluid drag is also disabled 167 y_top = -Inf 168 for grain in sim.grains 169 if y_top < grain.lin_pos[2] + grain.contact_radius 170 y_top = grain.lin_pos[2] + grain.contact_radius 171 end 172 end 173 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top, 174 bc="normal stress", 175 normal_stress=-20e3) 176 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass 177 @info("Placing top wall at y=$y_top") 178 179 # Resize the grid to span the current state 180 Granular.fitGridToGrains!(sim, sim.ocean) 181 182 Granular.setTimeStep!(sim) 183 Granular.setTotalTime!(sim, t_total) 184 Granular.setOutputFileInterval!(sim, file_dt) 185 render && Granular.plotGrains(sim, verbose=true) 186 187 ############################################################################ 188 #### RUN THE SIMULATION 189 ############################################################################ 190 191 Granular.run!(sim) 192 Granular.writeSimulation(sim) 193 194 render && Granular.render(sim, trim=false, animation=true) 195 end 196 197 function main() 198 parsed_args = parse_command_line() 199 report_args(parsed_args) 200 201 seed = parsed_args["seed"] 202 run_simulation(parsed_args["id"] * "-seed" * string(seed), 203 parsed_args["E"], 204 parsed_args["nu"], 205 parsed_args["mu_d"], 206 parsed_args["tensile_strength"], 207 parsed_args["tensile_strength_std_dev"], 208 parsed_args["shear_strength"], 209 parsed_args["fracture_toughness"], 210 parsed_args["r_min"], 211 parsed_args["r_max"], 212 parsed_args["thickness"], 213 parsed_args["rotating"], 214 parsed_args["compressive_velocity"], 215 seed) 216 end 217 218 main() 219 end