pipeline_test.jl (8164B)
1 import Granular 2 import JLD2 3 import PyPlot 4 import Dates 5 6 include("indent.jl") 7 8 #call this script as "@elapsed include("initialization.jl")" in order to time it 9 10 t_start = Dates.now() 11 12 t_init = 0.5 # duration of initialization [s] 13 t_stack = 0.4 #duration for each stacking 14 15 g = [0.,-9.8] 16 17 nx = 25 #125 18 ny = 8 #80 19 stacks = 0 20 21 ngrains = nx*ny*(stacks+1) 22 23 suffix = "_test" 24 25 r_min = 0.05 26 r_max = r_min*sqrt(2) #max grain size is double area of smallest grain size 27 28 SimSettings = Dict() 29 30 SimSettings["nx"] = nx 31 SimSettings["ny"] = ny 32 SimSettings["r_min"] = r_min 33 SimSettings["r_max"] = r_max 34 35 #JLD2.save("SimSettings$(ngrains)$(suffix).jld2", SimSettings) 36 37 gsd_type = "powerlaw" 38 gsd_powerlaw_exponent = -1.8 39 gsd_seed = 3 40 41 # mechanical properties of grains 42 youngs_modulus = 2e7 #elastic modulus 43 poissons_ratio = 0.185 #shear stiffness ratio 44 tensile_strength = 0.0 #strength of bonds between grains 45 contact_dynamic_friction = 0.4 #friction between grains 46 rotating = true #can grains rotate or not 47 48 sim = Granular.createSimulation(id="init") 49 sim.id = "init$(ngrains)$(suffix)" 50 51 Granular.regularPacking!(sim, #simulation object 52 [nx, ny], #number of grains along x and y axis 53 r_min, #smallest grain size 54 r_max, #largest grain size 55 tiling="triangular", #how the grains are tiled 56 origo = [0.0,0.0], 57 size_distribution=gsd_type, 58 size_distribution_parameter=gsd_powerlaw_exponent, 59 seed=gsd_seed) 60 61 for grain in sim.grains #go through all grains in sim 62 grain.youngs_modulus = youngs_modulus 63 grain.poissons_ratio = poissons_ratio 64 grain.tensile_strength = tensile_strength 65 grain.contact_dynamic_friction = contact_dynamic_friction 66 grain.rotating = rotating 67 end 68 69 Granular.fitGridToGrains!(sim, sim.ocean) 70 71 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 72 verbose=false) 73 #Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west") 74 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west") 75 76 77 for grain in sim.grains 78 Granular.addBodyForce!(grain, grain.mass*g) 79 Granular.disableOceanDrag!(grain) 80 grain.contact_viscosity_normal = 1e4 # N/(m/s) 81 end 82 83 Granular.setTimeStep!(sim) 84 85 Granular.setTotalTime!(sim, t_init) 86 87 Granular.setOutputFileInterval!(sim, .01) 88 89 Granular.run!(sim) 90 91 #Granular.writeSimulation(sim, 92 # filename = "init$(ngrains)$(suffix).jld2", 93 # folder = "init$(ngrains)$(suffix)") 94 95 #Stack it on top of each other 96 97 temp = deepcopy(sim) 98 99 for i = 1:stacks 100 101 global y_top = -Inf 102 for grain in sim.grains 103 if y_top < grain.lin_pos[2] #+ grain.contact_radius 104 global y_top = grain.lin_pos[2] #+ grain.contact_radius 105 end 106 end 107 108 109 for grain in temp.grains 110 111 lin_pos_lifted = [0.0,0.0] 112 113 lin_pos_lifted[1] = grain.lin_pos[1] 114 lin_pos_lifted[2] = grain.lin_pos[2] + y_top + r_max 115 116 117 Granular.addGrainCylindrical!(sim, 118 lin_pos_lifted, 119 grain.contact_radius, 120 grain.thickness, 121 youngs_modulus = grain.youngs_modulus, 122 poissons_ratio = grain.poissons_ratio, 123 tensile_strength = grain.tensile_strength, 124 contact_dynamic_friction = grain.contact_dynamic_friction, 125 verbose = false) 126 127 end 128 129 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) 130 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 131 verbose=false) 132 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east west", 133 verbose=false) 134 Granular.setTotalTime!(sim,t_stack) 135 Granular.setTimeStep!(sim) 136 Granular.setOutputFileInterval!(sim, .01) 137 138 for grain in sim.grains 139 Granular.addBodyForce!(grain, grain.mass*g) 140 Granular.disableOceanDrag!(grain) 141 end 142 143 #Granular.resetTime!(sim) 144 sim.time_iteration = 0 145 sim.time = 0.0 146 sim.file_time_since_output_file = 0. 147 148 Granular.setTotalTime!(sim, t_stack) 149 Granular.setTimeStep!(sim) 150 Granular.setOutputFileInterval!(sim, .01) 151 152 Granular.run!(sim) 153 154 end 155 156 # add a lower boundary consisting of grains bound together 157 # do this by creating a new simulation where the grains are bonded using 158 # findContacts! after creating overlapping grains. Then set their contact 159 # strengths to an infinite amount, but do not allow new contacts 160 161 carpet = Granular.createSimulation(id="init_carpet") 162 163 bot_r = r_min #radius of bottom layer grains 164 165 left_edge = round(sim.ocean.origo[1],digits=2) 166 length = round(sim.ocean.L[1],digits=2) 167 right_edge = left_edge+length 168 169 for i = left_edge+(bot_r/2):bot_r*1.99:left_edge+length 170 171 bot_pos = [i,round(sim.ocean.origo[2]-bot_r,digits=2)] 172 173 Granular.addGrainCylindrical!(carpet, 174 bot_pos, 175 bot_r, 176 0.1, 177 verbose = false, 178 tensile_strength = Inf, 179 shear_strength = Inf, 180 contact_stiffness_normal = Inf, 181 contact_stiffness_tangential = Inf) 182 183 end 184 185 Granular.findContactsAllToAll!(carpet) #find the grain contacts 186 187 #carpet.grains[10].fixed = true 188 #carpet.grains[10].lin_vel[1:2] = [0.0,0.0] 189 190 191 #add the carpet to the main simulation object 192 append!(sim.grains,carpet.grains) 193 194 195 #fit the ocean to the added grains and run to let the basin settle 196 197 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false) 198 199 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south", 200 verbose=false) 201 202 sim.time_iteration = 0 203 sim.time = 0.0 204 sim.file_time_since_output_file = 0. 205 Granular.setTotalTime!(sim, 0.2) 206 Granular.setTimeStep!(sim) 207 Granular.setOutputFileInterval!(sim, .01) 208 209 Granular.run!(sim) 210 211 Granular.writeSimulation(sim, 212 filename = "stacked$(ngrains)$(suffix).jld2") 213 214 ########## Add indenter ############# 215 216 temp_indent = createSimulation("id=temp_indent") 217 218 219 left_edge = round(sim.ocean.origo[1],digits=2) 220 length = round(sim.ocean.L[1],digits=2) 221 222 width = length/3 223 hw_ratio = 0.2 224 init_vertex_pos = [(length+left_edge)/2,-0.2] 225 grain_radius = 0.05 226 227 228 vertex_x = init_vertex_pos[1] 229 vertex_y = width*hw_ratio*sin((pi/width)*vertex_x) 230 231 232 for i = 0:grain_radius*2:width#manipulate the ocean grid 233 234 x_pos = i 235 236 y_pos = width*hw_ratio*sin(pi/width*x_pos) 237 238 Granular.addGrainCylindrical!(temp_indent, 239 [x_pos+init_vertex_pos[1]-width/2,y_pos+vertex_y+init_vertex_pos[2]], 240 grain_radius, 241 0.1, 242 fixed = true, 243 lin_vel = [0.0,0.5]) 244 end 245 246 append!(sim.grains,temp_indent.grains) 247 248 Granular.fitGridToGrains!(sim, 249 sim.ocean, 250 north_padding = 3.0, 251 verbose=false) 252 253 254 sim.time_iteration = 0 255 sim.time = 0.0 256 sim.file_time_since_output_file = 0. 257 Granular.setTotalTime!(sim, 0.5) 258 Granular.setTimeStep!(sim) 259 Granular.setOutputFileInterval!(sim, .01) 260 261 while sim.time < sim.time_total 262 for grain in carpet.grains 263 if grain.lin_vel[2] < 0 #&& grain.lin_pos[2] < r_min #find some lower boundary here? 264 grain.lin_vel[1:2] = [0.,0.] 265 end 266 end 267 Granular.run!(sim,single_step=true) 268 end 269 270 Granular.writeSimulation(sim, 271 filename = "indented$(ngrains)_test.jld2") 272 273 274 #print time elapsed 275 t_now = Dates.now() 276 dur = Dates.canonicalize(t_now-t_start) 277 print("Time elapsed: ",dur)