strait.jl (6879B)
1 #!/usr/bin/env julia 2 import Granular 3 using Random 4 5 sim = Granular.createSimulation(id="strait") 6 n = [10, 10, 2] 7 8 #sim = Granular.createSimulation(id="nares_strait_coarse_elast") 9 #n = [6, 6, 2] 10 11 # Initialize ocean 12 Lx = 50.e3 13 Lx_constriction = 5e3 14 L = [Lx, Lx*1.5, 1e3] 15 Ly_constriction = 20e3 16 sim.ocean = Granular.createRegularOceanGrid(n, L, name="poiseuille_flow") 17 sim.ocean.v[:, :, 1, 1] = 1e-8.*((sim.ocean.xq .- Lx/2.).^2.0 .- Lx^2.0/4.0) 18 19 # Initialize confining walls, which are grains that are fixed in space 20 r = minimum(L[1:2]/n[1:2])/2. 21 r_min = r/4. 22 h = 1. 23 24 ## N-S segments 25 r_walls = r_min 26 for y in range((L[2] - Ly_constriction)/2., 27 stop=Ly_constriction + (L[2] - Ly_constriction)/2.0, 28 length=Int(round(Ly_constriction/(r_walls*2)))) 29 Granular.addGrainCylindrical!(sim, [(Lx - Lx_constriction)/2.0, y], 30 r_walls, 31 h, fixed=true, verbose=false) 32 end 33 for y in range((L[2] - Ly_constriction)/2.0, 34 stop=Ly_constriction + (L[2] - Ly_constriction)/2.0, 35 length=Int(round(Ly_constriction/(r_walls*2)))) 36 Granular.addGrainCylindrical!(sim, 37 [Lx_constriction + 38 (L[1] - Lx_constriction)/2.0, y], r_walls, h, 39 fixed=true, verbose=false) 40 end 41 42 dx = 2.0*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction))) 43 44 ## NW diagonal 45 x = r_walls:dx:((Lx - Lx_constriction)/2.) 46 y = range(L[2] - r_walls, 47 stop=(L[2] - Ly_constriction)/2. + Ly_constriction + r_walls, 48 length=length(x)) 49 for i in 1:length(x) 50 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 51 verbose=false) 52 end 53 54 ## NE diagonal 55 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 56 y = range(L[2] - r_walls, 57 stop=(L[2] - Ly_constriction)/2. + Ly_constriction + r_walls, 58 length=length(x)) 59 for i in 1:length(x) 60 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 61 verbose=false) 62 end 63 64 ## SW diagonal 65 x = r_walls:dx:((Lx - Lx_constriction)/2.) 66 y = range(r, stop=(L[2] - Ly_constriction)/2.0 - r_walls, 67 length=length(x)) 68 for i in 1:length(x) 69 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 70 verbose=false) 71 end 72 73 ## SE diagonal 74 x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction) 75 y = range(r_walls, stop=(L[2] - Ly_constriction)/2. - r_walls, length=length(x)) 76 for i in 1:length(x) 77 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 78 verbose=false) 79 end 80 81 n_walls = length(sim.grains) 82 @info "added $(n_walls) fixed grains as walls" 83 84 # Initialize grains in wedge north of the constriction 85 dy = sqrt((2.0*r_walls)^2.0 - dx^2.0) 86 spacing_to_boundaries = 4.0*r 87 floe_padding = 0.5*r 88 noise_amplitude = floe_padding 89 Random.seed!(1) 90 let 91 iy = 1 92 for y in (L[2] - r - noise_amplitude):(-(2.0*r + floe_padding)):((L[2] - 93 Ly_constriction)/2.0 + Ly_constriction) 94 for x in (r + noise_amplitude):(2.0*r + floe_padding):(Lx - r - 95 noise_amplitude) 96 if iy % 2 == 0 97 x += 1.5*r 98 end 99 100 x_ = x + noise_amplitude*(0.5 - rand()) 101 y_ = y + noise_amplitude*(0.5 - rand()) 102 103 if y_ < -dy/dx*x_ + L[2] + spacing_to_boundaries 104 continue 105 end 106 107 if y_ < dy/dx*x_ + (L[2] - dy/dx*Lx) + spacing_to_boundaries 108 continue 109 end 110 111 r_rand = r_min + rand()*(r - r_min) 112 Granular.addGrainCylindrical!(sim, [x_, y_], r_rand, h, verbose=false) 113 end 114 iy += 1 115 end 116 end 117 n = length(sim.grains) - n_walls 118 @info "added $n grains" 119 120 # Remove old simulation files 121 Granular.removeSimulationFiles(sim) 122 123 k_n = 1e7 # N/m 124 k_t = k_n 125 #gamma_t = 1e7 # N/(m/s) 126 gamma_t = 0.0 127 mu_d = 0.7 128 rotating = true 129 for i=1:length(sim.grains) 130 sim.grains[i].contact_stiffness_normal = k_n 131 sim.grains[i].contact_stiffness_tangential = k_t 132 sim.grains[i].contact_viscosity_tangential = gamma_t 133 sim.grains[i].contact_dynamic_friction = mu_d 134 sim.grains[i].rotating = rotating 135 end 136 137 # Set temporal parameters 138 Granular.setTotalTime!(sim, 6.0*60.0*60.0) 139 Granular.setOutputFileInterval!(sim, 60.0) 140 Granular.setTimeStep!(sim) 141 142 # Run simulation for 10 time steps, then add new grains the top 143 while sim.time < sim.time_total 144 for it=1:10 145 Granular.run!(sim, status_interval=1, single_step=true) 146 end 147 for i=1:size(sim.ocean.xh, 1) 148 if sim.ocean.grain_list[i, end] == [] 149 150 x, y = Granular.getCellCenterCoordinates(sim.ocean.xh, 151 sim.ocean.yh, 152 i, size(sim.ocean.xh, 2)) 153 154 # Enable for high mass flux 155 r_rand = r_min + rand()*(r - r_min) 156 Granular.addGrainCylindrical!(sim, [x-r, y-r], r_rand, h, 157 verbose=false, 158 contact_stiffness_normal=k_n, 159 contact_stiffness_tangential=k_t, 160 contact_viscosity_tangential=gamma_t, 161 contact_dynamic_friction = mu_d, 162 rotating=rotating) 163 r_rand = r_min + rand()*(r - r_min) 164 Granular.addGrainCylindrical!(sim, [x+r, y-r], r_rand, h, 165 verbose=false, 166 contact_stiffness_normal=k_n, 167 contact_stiffness_tangential=k_t, 168 contact_viscosity_tangential=gamma_t, 169 contact_dynamic_friction = mu_d, 170 rotating=rotating) 171 r_rand = r_min + rand()*(r - r_min) 172 Granular.addGrainCylindrical!(sim, [x+r, y+r], r_rand, h, 173 verbose=false, 174 contact_stiffness_normal=k_n, 175 contact_stiffness_tangential=k_t, 176 contact_viscosity_tangential=gamma_t, 177 contact_dynamic_friction = mu_d, 178 rotating=rotating) 179 r_rand = r_min + rand()*(r - r_min) 180 Granular.addGrainCylindrical!(sim, [x-r, y+r], r_rand, h, 181 verbose=false, 182 contact_stiffness_normal=k_n, 183 contact_stiffness_tangential=k_t, 184 contact_viscosity_tangential=gamma_t, 185 contact_dynamic_friction = mu_d, 186 rotating=rotating) 187 188 # Enable for low mass flux 189 #x += noise_amplitude*(0.5 - rand()) 190 #y += noise_amplitude*(0.5 - rand()) 191 #Granular.addGrainCylindrical!(sim, [x, y], r, h, verbose=false) 192 end 193 end 194 end