layer_basin.jl (3778B)
1 import Granular 2 import JLD2 3 import PyPlot 4 import Dates 5 6 id = "simulation40000" # id of simulation to load, just write the folder 7 # name here 8 9 # Layer interface positions 10 # defined as decimals 11 # ex: [0.4,0.6,1] will give two boundaries at 0.4 from the bottom 12 # and 0.6 from the bottom. ie 3 layers 13 # include 1 in the end as roof. 14 interfaces = [0,0.4,0.6,1] 15 16 # mechanical properties for each layer 17 youngs_modulus = [2e7,2e7,2e7] # elastic modulus 18 poissons_ratio = [0.185,0.185,0.185] # shear stiffness ratio 19 tensile_strength = [0.3,0.01,0.3] # strength of bonds between grains 20 shear_strength = [0.3,0.01,0.3] # shear stregth of bonds 21 contact_dynamic_friction = [0.4,0.01,0.4] # friction between grains 22 rotating = [true,true,true] # can grains rotate or not 23 color = [1,2,1] 24 25 #carpet_youngs_modulus = 2e7 26 #carpet_poissons_ratio = 0.185 27 #carpet_tensile_strength = 1e16 28 #carpet_contact_dynamic_friction = 0.4 29 #carpet_rotating = true 30 #carpet_shear_strength = 1e16 31 32 sim = Granular.readSimulation("$(id)/comp.jld2") 33 SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2") 34 35 sim.walls = Granular.WallLinearFrictionless[] # remove existing walls 36 37 Granular.zeroKinematics!(sim) # end any movement 38 39 y_top = -Inf 40 for grain in sim.grains 41 grain.contact_viscosity_normal = 0 42 if y_top < grain.lin_pos[2] + grain.contact_radius 43 global y_top = grain.lin_pos[2] + grain.contact_radius 44 end 45 end 46 47 y_bot = Inf 48 for grain in sim.grains 49 if y_bot > grain.lin_pos[2] - grain.contact_radius 50 global y_bot = grain.lin_pos[2] - grain.contact_radius 51 end 52 end 53 54 # quick fix to make the color = 0 flag for grains belonging to the carpet. 55 # this should be done in the newer versions of init_basin.jl instead 56 57 for grain in sim.grains 58 if grain.lin_pos[2] == -0.05 59 grain.color = 0 60 else 61 grain.color = 1 62 end 63 end 64 65 66 h = y_top-y_bot #depth of basin 67 68 interfaces *= h 69 70 for grain in sim.grains 71 72 for i = 2:size(interfaces,1) 73 74 if grain.lin_pos[2] <= interfaces[i] && grain.lin_pos[2] > interfaces[i-1] && grain.color != 0 75 76 grain.youngs_modulus = youngs_modulus[i-1] 77 grain.poissons_ratio = poissons_ratio[i-1] 78 grain.tensile_strength = tensile_strength[i-1] 79 grain.shear_strength = shear_strength[i-1] 80 grain.contact_dynamic_friction = contact_dynamic_friction[i-1] 81 grain.rotating = rotating[i-1] 82 grain.color = color[i-1] 83 84 end 85 end 86 end 87 88 # Create the bonds between grains by expanding all grains by a small amount 89 # then search and establish contacts and then reduce the size of the grains again 90 91 size_increasing_factor = 1.10 # factor by which contact radius should be increased 92 # to search for contacts 93 increase_array = [] 94 95 #increase the contact radius 96 for grain in sim.grains 97 if grain.color != 0 98 contact_radius_increase = (grain.contact_radius*size_increasing_factor)-grain.contact_radius 99 grain.contact_radius += contact_radius_increase 100 append!(increase_array,contact_radius_increase) 101 elseif grain.color == 0 102 append!(increase_array,0) 103 end 104 end 105 106 Granular.findContacts!(sim,method="ocean grid") 107 #Granular.findContactsAllToAll!(sim) # find the grain contacts 108 #Granular.run!(sim,single_step=true) 109 110 #reduce the contact radius again 111 for i = 1:size(sim.grains,1) 112 sim.grains[i].contact_radius -= increase_array[i] 113 end 114 115 116 cd("$id") 117 sim.id = "layered" 118 119 #Granular.resetTime!(sim) 120 #Granular.setTotalTime!(sim,0.5) 121 #Granular.run!(sim) 122 123 Granular.resetTime!(sim) 124 Granular.run!(sim,single_step=true) 125 126 cd("..") 127 128 Granular.writeSimulation(sim, 129 filename = "$(id)/layered.jld2")