tshear.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
---
tshear.jl (10814B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import JLD2
5 import PyPlot
6
7 ################################################################################
8 #### SIMULATION PARAMETERS #
9 ################################################################################
10 let
11
12 # Common simulation identifier
13 id_prefix = "test0"
14
15 # Gravitational acceleration vector (cannot be zero; required for Step 1)
16 g = [0., -9.8]
17
18 # Grain package geometry during initialization
19 nx = 10 # Grains along x (horizontal)
20 ny = 50 # Grains along y (vertical)
21
22 # Grain-size parameters
23 r_min = 0.03 # Min. grain radius [m]
24 r_max = 0.1 # Max. grain radius [m]
25 gsd_type = "powerlaw" # "powerlaw" or "uniform" sizes between r_min and r_max
26 gsd_powerlaw_exponent = -1.8 # GSD power-law exponent
27 gsd_seed = 1 # Value to seed random-size generation
28
29 # Grain mechanical properties
30 youngs_modulus = 2e7 # Elastic modulus [Pa]
31 poissons_ratio = 0.185 # Shear-stiffness ratio [-]
32 tensile_strength = 0.0 # Inter-grain bond strength [Pa]
33 contact_dynamic_friction = 0.4 # Coulomb-frictional coefficient [-]
34 rotating = true # Allow grain rotation
35
36 # Normal stress for the consolidation and shear [Pa]
37 N = 20e3
38
39 # Shear velocity to apply to the top grains [m/s]
40 vel_shear = 0.5
41
42 # Simulation duration of individual steps [s]
43 t_init = 2.0
44 t_cons = 2.5
45 t_shear = 5.0
46
47 ################################################################################
48 #### Step 1: Create a loose granular assemblage and let it settle at -y #
49 ################################################################################
50 sim = Granular.createSimulation(id="$(id_prefix)-init")
51
52 Granular.regularPacking!(sim, [nx, ny], r_min, r_max,
53 size_distribution=gsd_type,
54 size_distribution_parameter=gsd_powerlaw_exponent,
55 seed=gsd_seed)
56
57 # Set grain mechanical properties
58 for grain in sim.grains
59 grain.youngs_modulus = youngs_modulus
60 grain.poissons_ratio = poissons_ratio
61 grain.tensile_strength = tensile_strength
62 grain.contact_dynamic_friction = contact_dynamic_friction
63 grain.rotating = rotating
64 end
65
66 # Create a grid for contact searching spanning the extent of the grains
67 Granular.fitGridToGrains!(sim, sim.ocean)
68
69 # Make the top and bottom boundaries impermeable, and the side boundaries
70 # periodic, which will come in handy during shear
71 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "north south",
72 verbose=false)
73 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "east west")
74
75 # Add gravitational acceleration to all grains and disable ocean-grid drag.
76 # Also add viscous energy dissipation between grains, which is disabled before
77 # consolidation and shear.
78 for grain in sim.grains
79 Granular.addBodyForce!(grain, grain.mass*g)
80 Granular.disableOceanDrag!(grain)
81 grain.contact_viscosity_normal = 1e4 # N/(m/s)
82 end
83
84 # Automatically set the computational time step based on grain sizes and
85 # properties
86 Granular.setTimeStep!(sim)
87
88 # Set the total simulation time for this step [s]
89 # This value may need tweaking if grain sizes or numbers are adjusted.
90 Granular.setTotalTime!(sim, t_init)
91
92 # Set the interval in model time between simulation files [s]
93 Granular.setOutputFileInterval!(sim, .02)
94
95 # Visualize the grain-size distribution
96 Granular.plotGrainSizeDistribution(sim)
97
98 # Start the simulation
99 Granular.run!(sim)
100
101 # Try to render the simulation if `pvpython` is installed on the system
102 Granular.render(sim, trim=false)
103
104 # Save the simulation state to disk in case we need to reuse the current state
105 Granular.writeSimulation(sim)
106
107 # Also copy the simulation in memory, in case we want to loop over different
108 # normal stresses below:
109 sim_init = deepcopy(sim)
110
111
112 ################################################################################
113 #### Step 2: Consolidate the previous output under a constant normal stress #
114 ################################################################################
115
116 # Rename the simulation so we don't overwrite output from the previous step
117 sim.id = "$(id_prefix)-cons-N$(N)Pa"
118
119 # Set all linear and rotational velocities to zero
120 Granular.zeroKinematics!(sim)
121
122 # Add a dynamic wall to the top which adds a normal stress downwards. The
123 # normal of this wall is downwards, and we place it at the top of the granular
124 # assemblage. Here, the inter-grain viscosity is also removed.
125 y_top = -Inf
126 for grain in sim.grains
127 grain.contact_viscosity_normal = 0.
128 if y_top < grain.lin_pos[2] + grain.contact_radius
129 y_top = grain.lin_pos[2] + grain.contact_radius
130 end
131 end
132 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
133 bc="normal stress", normal_stress=-N,
134 contact_viscosity_normal=1e3)
135 @info "Placing top wall at y=$y_top"
136
137 # Resize the grid to span the current state
138 Granular.fitGridToGrains!(sim, sim.ocean)
139
140 # Lock the grains at the very bottom so that the lower boundary is rough
141 y_bot = Inf
142 for grain in sim.grains
143 if y_bot > grain.lin_pos[2] - grain.contact_radius
144 y_bot = grain.lin_pos[2] - grain.contact_radius
145 end
146 end
147 fixed_thickness = 2. * r_max
148 for grain in sim.grains
149 if grain.lin_pos[2] <= fixed_thickness
150 grain.fixed = true # set x and y acceleration to zero
151 end
152 end
153
154 # Set current time to zero and reset output file counter
155 Granular.resetTime!(sim)
156
157 # Set the simulation time to run the consolidation for
158 Granular.setTotalTime!(sim, t_cons)
159
160 # Run the consolidation experiment, and monitor top wall position over time
161 time = Float64[]
162 compaction = Float64[]
163 effective_normal_stress = Float64[]
164 while sim.time < sim.time_total
165
166 for i=1:100 # run for 100 steps before measuring shear stress and dilation
167 Granular.run!(sim, single_step=true)
168 end
169
170 append!(time, sim.time)
171 append!(compaction, sim.walls[1].pos)
172 append!(effective_normal_stress, Granular.getWallNormalStress(sim))
173
174 end
175 defined_normal_stress = ones(length(effective_normal_stress)) *
176 Granular.getWallNormalStress(sim, stress_type="effective")
177 PyPlot.subplot(211)
178 PyPlot.subplots_adjust(hspace=0.0)
179 ax1 = PyPlot.gca()
180 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels
181 PyPlot.plot(time, compaction)
182 PyPlot.ylabel("Top wall height [m]")
183 PyPlot.subplot(212, sharex=ax1)
184 PyPlot.plot(time, defined_normal_stress)
185 PyPlot.plot(time, effective_normal_stress)
186 PyPlot.xlabel("Time [s]")
187 PyPlot.ylabel("Normal stress [Pa]")
188 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
189 PyPlot.clf()
190
191 # Try to render the simulation if `pvpython` is installed on the system
192 Granular.render(sim, trim=false)
193
194 # Save the simulation state to disk in case we need to reuse the consolidated
195 # state (e.g. different shear velocities below)
196 Granular.writeSimulation(sim)
197
198 # Also copy the simulation in memory, in case we want to loop over different
199 # normal stresses below:
200 sim_cons = deepcopy(sim)
201
202
203 ################################################################################
204 #### Step 3: Shear the consolidated assemblage with a constant velocity #
205 ################################################################################
206
207 # Rename the simulation so we don't overwrite output from the previous step
208 sim.id = "$(id_prefix)-shear-N$(N)Pa-vel_shear$(vel_shear)m-s"
209
210 # Set all linear and rotational velocities to zero
211 Granular.zeroKinematics!(sim)
212
213 # Set current time to zero and reset output file counter
214 Granular.resetTime!(sim)
215
216 # Set the simulation time to run the shear experiment for
217 Granular.setTotalTime!(sim, t_shear)
218
219 # Run the shear experiment
220 time = Float64[]
221 shear_stress = Float64[]
222 shear_strain = Float64[]
223 dilation = Float64[]
224 thickness_initial = sim.walls[1].pos - y_bot
225 x_min = +Inf
226 x_max = -Inf
227 for grain in sim.grains
228 if x_min > grain.lin_pos[1] - grain.contact_radius
229 x_min = grain.lin_pos[1] - grain.contact_radius
230 end
231 if x_max < grain.lin_pos[1] + grain.contact_radius
232 x_max = grain.lin_pos[1] + grain.contact_radius
233 end
234 end
235 surface_area = (x_max - x_min)
236 shear_force = 0.
237 while sim.time < sim.time_total
238
239 # Prescribe the shear velocity to the uppermost grains
240 for grain in sim.grains
241 if grain.lin_pos[2] >= sim.walls[1].pos - fixed_thickness
242 grain.fixed = true
243 grain.allow_y_acc = true
244 grain.lin_vel[1] = vel_shear
245 elseif grain.lin_pos[2] > fixed_thickness # do not free bottom grains
246 grain.fixed = false
247 end
248 end
249
250 for i=1:100 # run for 100 steps before measuring shear stress and dilation
251 Granular.run!(sim, single_step=true)
252 end
253
254 append!(time, sim.time)
255
256 # Determine the current shear stress
257 shear_force = 0.
258 for grain in sim.grains
259 if grain.fixed && grain.allow_y_acc
260 shear_force += -grain.force[1]
261 end
262 end
263 append!(shear_stress, shear_force/surface_area)
264
265 # Determine the current shear strain
266 append!(shear_strain, sim.time*vel_shear/thickness_initial)
267
268 # Determine the current dilation
269 append!(dilation, (sim.walls[1].pos - y_bot)/thickness_initial)
270
271 end
272
273 # Try to render the simulation if `pvpython` is installed on the system
274 Granular.render(sim, trim=false)
275
276 # Save the simulation state to disk in case we need to reuse the sheared state
277 Granular.writeSimulation(sim)
278
279 # Plot time vs. shear stress and dilation
280 PyPlot.subplot(211)
281 PyPlot.subplots_adjust(hspace=0.0)
282 ax1 = PyPlot.gca()
283 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels
284 PyPlot.plot(time, shear_stress)
285 PyPlot.ylabel("Shear stress [Pa]")
286 PyPlot.subplot(212, sharex=ax1)
287 PyPlot.plot(time, dilation)
288 PyPlot.xlabel("Time [s]")
289 PyPlot.ylabel("Volumetric strain [-]")
290 PyPlot.savefig(sim.id * "-time_vs_shear-stress.pdf")
291 PyPlot.clf()
292
293 # Plot shear strain vs. shear stress and dilation
294 PyPlot.subplot(211)
295 PyPlot.subplots_adjust(hspace=0.0)
296 ax1 = PyPlot.gca()
297 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels
298 PyPlot.plot(shear_strain, shear_stress)
299 PyPlot.ylabel("Shear stress [Pa]")
300 PyPlot.subplot(212, sharex=ax1)
301 PyPlot.plot(shear_strain, dilation)
302 PyPlot.xlabel("Shear strain [-]")
303 PyPlot.ylabel("Volumetric strain [-]")
304 PyPlot.savefig(sim.id * "-shear-strain_vs_shear-stress.pdf")
305 PyPlot.clf()
306
307 # Plot time vs. shear strain (boring when the shear experiment is rate controlled)
308 PyPlot.plot(time, shear_strain)
309 PyPlot.xlabel("Time [s]")
310 PyPlot.ylabel("Shear strain [-]")
311 PyPlot.savefig(sim.id * "-time_vs_shear-strain.pdf")
312 PyPlot.clf()
313 end