Introduction
Introduction Statistics Contact Development Disclaimer Help
tpipeline_test.jl - granular-basin - tectonic deformation experiments with Gran…
git clone git://src.adamsgaard.dk/granular-basin
Log
Files
Refs
README
---
tpipeline_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…
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 s…
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_exponen…
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 sou…
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_modulu…
122 poissons_ratio = grain.poissons_rati…
123 tensile_strength = grain.tensile_str…
124 contact_dynamic_friction = grain.con…
125 verbose = false)
126
127 end
128
129 Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)
130 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north…
131 verb…
132 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "east …
133 verb…
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 sou…
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_…
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 …
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)
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.