init_basin.jl (9058B)
1 import Granular 2 import JLD2 3 import PyPlot 4 import Dates 5 6 t_start = Dates.now() # Save the start time, print the end time later. 7 8 ############# Initialization Settings ############# 9 10 t_init = 1.5 # duration of initialization [s] 11 t_stack = 1.0 # duration for each stack to settle [s] 12 13 g = [0.,-9.8] # vector for direction and magnitude of gravitational acceleration of grains 14 15 ngrains = 250 # total number of grains 16 aspect_ratio = 4 # should be x times as wide as it is tall 17 18 mkpath("simulation$(ngrains)") 19 20 21 stacks = 2 # number of duplicate stacks on top of the initial grains 22 23 ny = sqrt((ngrains)/aspect_ratio) # number of grain rows 24 nx = aspect_ratio*ny*(stacks+1) # number of grain columns 25 26 ny = Int(round(ny/(stacks+1))) 27 nx = Int(round(nx/(stacks+1))) 28 29 30 r_min = 0.05 # minimum radius of grains 31 r_max = r_min*sqrt(2) # max radius of grains, double the area of minimum 32 33 # Grain-size distribution parameters 34 gsd_type = "powerlaw" # type grain-size distribution 35 gsd_powerlaw_exponent = -1.8 # exponent if powerlaw is used for grain-size distribution 36 gsd_seed = 3 # random seed for the distribution of grain-sizes 37 38 # Initial mechanical properties of grains 39 youngs_modulus = 2e7 # elastic modulus 40 poissons_ratio = 0.185 # shear stiffness ratio 41 tensile_strength = 0.0 # strength of bonds between grains 42 contact_dynamic_friction = 0.4 # friction between grains 43 rotating = true # can grains rotate or not 44 45 # Save some settings in a dictionary 46 SimSettings = Dict() 47 SimSettings["nx"] = nx 48 SimSettings["ny"] = ny 49 SimSettings["stacks"] = stacks 50 SimSettings["ngrains"] = ngrains 51 SimSettings["r_min"] = r_min 52 SimSettings["r_max"] = r_max 53 54 # this section has been moved further up 55 56 ############# Initialize simulation and grains ############# 57 58 sim = Granular.createSimulation(id="simulation$(ngrains)") 59 60 Granular.regularPacking!(sim, #simulation object 61 [nx, ny], #number of grains along x and y axis 62 r_min, #smallest grain size 63 r_max, #largest grain size 64 tiling="triangular", #how the grains are tiled 65 origo = [0.0,0.0], 66 size_distribution=gsd_type, 67 size_distribution_parameter=gsd_powerlaw_exponent, 68 seed=gsd_seed, 69 color = 1) 70 71 72 # set the indicated mechanical parameters for all grains 73 for grain in sim.grains 74 grain.youngs_modulus = youngs_modulus 75 grain.poissons_ratio = poissons_ratio 76 grain.tensile_strength = tensile_strength 77 grain.contact_dynamic_friction = contact_dynamic_friction 78 grain.rotating = rotating 79 end 80 81 # fit the ocean grid to the grains 82 Granular.fitGridToGrains!(sim, sim.ocean) 83 84 # set the boundary conditions 85 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 86 verbose=false) 87 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west", 88 verbose=false) 89 90 # add gravity to the grains and remove ocean drag 91 for grain in sim.grains 92 Granular.addBodyForce!(grain, grain.mass*g) 93 Granular.disableOceanDrag!(grain) 94 grain.contact_viscosity_normal = 1e4 # N/(m/s) 95 end 96 97 # run the simulation 98 Granular.setTimeStep!(sim) # set appropriate time steps 99 Granular.setTotalTime!(sim, t_init) # set total time 100 Granular.setOutputFileInterval!(sim, .01) # how often vtu files should be outputted 101 102 cd("simulation$(ngrains)") 103 104 sim.id = "init" 105 106 Granular.run!(sim) 107 108 109 110 ############# Stack the initialized grains ############# 111 112 temp = deepcopy(sim) 113 114 for i = 1:stacks 115 116 # find the position of the uppermost grain 117 global y_top = -Inf 118 119 for grain in sim.grains 120 if y_top < grain.lin_pos[2] 121 global y_top = grain.lin_pos[2] 122 end 123 end 124 125 # add duplicate grains above the initialized grains 126 for grain in temp.grains 127 128 lin_pos_lifted = [0.0,0.0] # preallocation of position of a 'lifted' grain 129 lin_pos_lifted[1] = grain.lin_pos[1] # x position of duplicate grain 130 lin_pos_lifted[2] = grain.lin_pos[2] + y_top + r_max # y-position of duplicate grain 131 132 133 Granular.addGrainCylindrical!(sim, 134 lin_pos_lifted, 135 grain.contact_radius, 136 grain.thickness, 137 youngs_modulus = grain.youngs_modulus, 138 poissons_ratio = grain.poissons_ratio, 139 tensile_strength = grain.tensile_strength, 140 contact_dynamic_friction = grain.contact_dynamic_friction, 141 verbose = false) 142 143 end 144 145 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) 146 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 147 verbose=false) 148 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west", 149 verbose=false) 150 Granular.setTotalTime!(sim,t_stack) 151 Granular.setTimeStep!(sim) 152 Granular.setOutputFileInterval!(sim, .01) 153 154 for grain in sim.grains 155 Granular.addBodyForce!(grain, grain.mass*g) 156 Granular.disableOceanDrag!(grain) 157 end 158 159 160 # instead of Granular.resetTime!(sim), the time parameters are reset manually, in order 161 # to allow for the continuation of vtu-file index from earlier 162 sim.time_iteration = 0 163 sim.time = 0.0 164 sim.file_time_since_output_file = 0. 165 166 Granular.setTotalTime!(sim, t_stack) 167 Granular.setTimeStep!(sim) 168 Granular.setOutputFileInterval!(sim, .01) 169 170 Granular.run!(sim) # let the duplicated grains settle. 171 172 end 173 174 175 176 177 ############# Lay a carpet ############# 178 179 carpet = Granular.createSimulation(id="init_carpet") # new simulation object for the carpet 180 181 bot_r = r_min # radius of carpet grains 182 183 left_edge = round(sim.ocean.origo[1],digits=2) # west edge of the carpet 184 length = round(sim.ocean.L[1],digits=2) # width of the carpet 185 right_edge = left_edge+length # east edge of the carpet 186 187 188 # Now loop over the carpet grain positions, the loop will create grains that overlap slightly 189 # in order to create the bonds needed 190 # color = 0 is used as a flag for the grains in the carpet 191 for i = left_edge+(bot_r/2):bot_r*1.999:left_edge+length 192 193 bot_pos = [i,round(sim.ocean.origo[2]-bot_r,digits=2)] # position of grain 194 195 Granular.addGrainCylindrical!(carpet, 196 bot_pos, 197 bot_r, 198 0.1, 199 verbose = false, 200 tensile_strength = Inf, 201 shear_strength = Inf, 202 #contact_stiffness_normal = Inf, 203 #contact_stiffness_tangential = Inf, 204 fixed = false, 205 color = 0) 206 end 207 208 #Granular.fitGridToGrains!(carpet,carpet.ocean,verbose=false) 209 210 211 212 Granular.findContactsAllToAll!(carpet) # find the grain contacts 213 214 215 216 append!(sim.grains,carpet.grains) # add the carpet grains to the main simulation object 217 # since the assignment will point to the carpet object, changes made to the carpet 218 # object will appear in the main simulation object 219 220 221 """ 222 #reset the grain contacts and make them very old 223 224 for grain in sim.grains 225 grain.contacts[:] .= 0 226 grain.n_contacts = 0 227 end 228 229 230 for grain in sim.grains 231 for ic=1:size(grain.contact_age,1) 232 grain.contact_age[ic] = 1e16 233 end 234 grain.strength_heal_rate = 1 # new bond stengthening 235 end 236 """ 237 238 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) # fit the ocean to the added grains 239 240 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 241 verbose=false) 242 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west", 243 verbose=false) 244 245 #Granular.findContacts!(sim,method="ocean grid") 246 247 # run the simulation shortly, to let the stacked grains settle on the carpet 248 sim.time_iteration = 0 249 sim.time = 0.0 250 sim.file_time_since_output_file = 0. 251 Granular.setTotalTime!(sim, 0.5) 252 Granular.setTimeStep!(sim) 253 Granular.setOutputFileInterval!(sim, .01) 254 255 256 Granular.run!(sim) 257 258 259 # save the simulation and the carpet objects 260 261 cd("..") 262 263 Granular.writeSimulation(sim, 264 filename = "simulation$(ngrains)/init.jld2") 265 266 JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings) 267 268 # print time elapsed 269 t_now = Dates.now() 270 dur = Dates.canonicalize(t_now-t_start) 271 print("Time elapsed: ",dur)