Introduction
Introduction Statistics Contact Development Disclaimer Help
tstrait.jl - Granular.jl - Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log
Files
Refs
README
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/…
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_w…
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=t…
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=t…
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=t…
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=le…
76 for i in 1:length(x)
77 Granular.addGrainCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=t…
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_amplitud…
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=…
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.x…
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=fa…
192 end
193 end
194 end
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.