tstrait.jl - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tstrait.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