tsimulation_shear.jl - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tsimulation_shear.jl (8959B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import JLD
5 import PyPlot
6 import ArgParse
7
8 const render = false
9
10 function parse_command_line()
11 s = ArgParse.ArgParseSettings()
12 ArgParse.@add_arg_table s begin
13 "--E"
14 help = "Young's modulus [Pa]"
15 arg_type = Float64
16 default = 2e7
17 "--nu"
18 help = "Poisson's ratio [-]"
19 arg_type = Float64
20 default = 0.285
21 "--mu_d"
22 help = "dynamic friction coefficient [-]"
23 arg_type = Float64
24 default = 0.3
25 "--tensile_strength"
26 help = "the maximum tensile strength [Pa]"
27 arg_type = Float64
28 default = 0.0
29 "--r_min"
30 help = "minimum grain radius [m]"
31 arg_type = Float64
32 default = 5.0
33 "--r_max"
34 help = "maximum grain radius [m]"
35 arg_type = Float64
36 default = 1.35e3
37 default = 50.0
38 "--thickness"
39 help = "grain thickness [m]"
40 arg_type = Float64
41 default = 1.0
42 "--rotating"
43 help = "allow the grains to rotate"
44 arg_type = Bool
45 default = true
46 "--shearvel"
47 help = "shear velocity [m/s]"
48 arg_type = Float64
49 default = 1.0
50 "--seed"
51 help = "seed for the pseudo RNG"
52 arg_type = Int
53 default = 1
54 "id"
55 help = "simulation id"
56 required = true
57 end
58 return ArgParse.parse_args(s)
59 end
60
61 function report_args(parsed_args)
62 println("Parsed args:")
63 for (arg,val) in parsed_args
64 println(" $arg => $val")
65 end
66 end
67
68 function run_simulation(id::String,
69 E::Float64,
70 nu::Float64,
71 mu_d::Float64,
72 tensile_strength::Float64,
73 r_min::Float64,
74 r_max::Float64,
75 thickness::Float64,
76 rotating::Bool,
77 shearvel::Float64,
78 seed::Int)
79
80 ############################################################################
81 #### SIMULATION PARAMETERS #
82 ############################################################################
83
84 # Common simulation identifier
85 const id_prefix = id
86
87 # Shear velocity to apply to the top grains [m/s]
88 const vel_shear = shearvel
89
90 # Normal stresses for consolidation and shear [Pa]
91 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3]
92
93 # Simulation duration of individual steps [s]
94 const t_shear = 800.0
95
96 tau_p = Float64[]
97 tau_u = Float64[]
98
99 for N in N_list
100
101 ########################################################################
102 #### Step 3: Shear the consolidated assemblage with a constant velocity#
103 ########################################################################
104
105 sim_cons = Granular.createSimulation("$(id_prefix)-cons-N$(N)Pa")
106 sim = Granular.readSimulation(sim_cons)
107
108 sim.id = "$(id_prefix)-shear-N$(N)Pa-vel_shear$(vel_shear)m-s"
109 Granular.removeSimulationFiles(sim)
110
111 # Set all linear and rotational velocities to zero
112 Granular.zeroKinematics!(sim)
113
114 # Set current time to zero and reset output file counter
115 Granular.resetTime!(sim)
116
117 # Set the simulation time to run the shear experiment for
118 Granular.setTotalTime!(sim, t_shear)
119
120 y_bot = Inf
121 for grain in sim.grains
122 if y_bot > grain.lin_pos[2] - grain.contact_radius
123 y_bot = grain.lin_pos[2] - grain.contact_radius
124 end
125 Granular.disableOceanDrag!(grain)
126 end
127 # Run the shear experiment
128 time = Float64[]
129 shear_stress = Float64[]
130 shear_strain = Float64[]
131 dilation = Float64[]
132 const thickness_initial = sim.walls[1].pos - y_bot
133 x_min = +Inf
134 x_max = -Inf
135 for grain in sim.grains
136 if x_min > grain.lin_pos[1] - grain.contact_radius
137 x_min = grain.lin_pos[1] - grain.contact_radius
138 end
139 if x_max < grain.lin_pos[1] + grain.contact_radius
140 x_max = grain.lin_pos[1] + grain.contact_radius
141 end
142 end
143 const fixed_thickness = 2.0*r_max
144 const surface_area = (x_max - x_min)
145 shear_force = 0.
146 while sim.time < sim.time_total
147
148 # Prescribe the shear velocity to the uppermost grains
149 for grain in sim.grains
150 if grain.lin_pos[2] >= sim.walls[1].pos - fixed_thickness
151 grain.fixed = true
152 grain.color = 1
153 grain.allow_y_acc = true
154 grain.lin_vel[1] = vel_shear
155 elseif grain.lin_pos[2] > fixed_thickness # do not free bottom
156 grain.fixed = false
157 grain.color = 0
158 end
159 end
160
161 for i=1:100
162 Granular.run!(sim, single_step=true)
163 end
164
165 append!(time, sim.time)
166
167 # Determine the current shear stress
168 shear_force = 0.0
169 for grain in sim.grains
170 if grain.fixed && grain.allow_y_acc
171 shear_force += -grain.force[1]
172 end
173 end
174 append!(shear_stress, shear_force/surface_area)
175
176 # Determine the current shear strain
177 append!(shear_strain, sim.time*vel_shear/thickness_initial)
178
179 # Determine the current dilation
180 append!(dilation, (sim.walls[1].pos - y_bot)/thickness_initial)
181
182 end
183
184 # Try to render the simulation if `pvpython` is installed on the system
185 if render
186 Granular.render(sim, trim=false)
187 end
188
189 #Granular.writeSimulation(sim)
190
191 # Plot time vs. shear stress and dilation
192 PyPlot.subplot(211)
193 PyPlot.subplots_adjust(hspace=0.0)
194 ax1 = PyPlot.gca()
195 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
196 PyPlot.plot(time, shear_stress)
197 PyPlot.ylabel("Shear stress [Pa]")
198 PyPlot.subplot(212, sharex=ax1)
199 PyPlot.plot(time, dilation)
200 PyPlot.xlabel("Time [s]")
201 PyPlot.ylabel("Volumetric strain [-]")
202 PyPlot.savefig(sim.id * "-time_vs_shear-stress.pdf")
203 PyPlot.clf()
204
205 # Plot shear strain vs. shear stress and dilation
206 PyPlot.subplot(211)
207 PyPlot.subplots_adjust(hspace=0.0)
208 ax1 = PyPlot.gca()
209 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
210 PyPlot.plot(shear_strain, shear_stress)
211 PyPlot.ylabel("Shear stress [Pa]")
212 PyPlot.subplot(212, sharex=ax1)
213 PyPlot.plot(shear_strain, dilation)
214 PyPlot.xlabel("Shear strain [-]")
215 PyPlot.ylabel("Volumetric strain [-]")
216 PyPlot.savefig(sim.id * "-shear-strain_vs_shear-stress.pdf")
217 PyPlot.clf()
218
219 # Plot time vs. shear strain (boring when the shear experiment is rate
220 # controlled)
221 PyPlot.plot(time, shear_strain)
222 PyPlot.xlabel("Time [s]")
223 PyPlot.ylabel("Shear strain [-]")
224 PyPlot.savefig(sim.id * "-time_vs_shear-strain.pdf")
225 PyPlot.clf()
226
227 writedlm(sim.id * "-data.txt", [time, shear_strain,
228 shear_stress, dilation])
229
230 append!(tau_p, maximum(shear_stress[1:100]))
231 append!(tau_u, mean(shear_stress[end-50:end]))
232 end
233
234 # Plot normal stress vs. peak/ultimate shear stress
235 writedlm(id_prefix * "-data.txt", [N_list, tau_p, tau_u])
236 PyPlot.subplot(211)
237 PyPlot.subplots_adjust(hspace=0.0)
238 ax1 = PyPlot.gca()
239 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x tick labels
240 PyPlot.plot(N_list, tau_p)
241 PyPlot.ylabel("Peak shear stress [Pa]")
242 PyPlot.subplot(212, sharex=ax1)
243 PyPlot.plot(N_list, tau_u)
244 PyPlot.xlabel("Shear strain [-]")
245 PyPlot.ylabel("Ultimate shear stress [Pa]")
246 PyPlot.xlim(0.-, maximum(N_list)*1.1)
247 PyPlot.savefig(id_prefix * "_mc.pdf")
248 PyPlot.clf()
249 end
250
251 function main()
252 parsed_args = parse_command_line()
253 report_args(parsed_args)
254
255 seed = parsed_args["seed"]
256 run_simulation(parsed_args["id"] * "-seed" * string(seed),
257 parsed_args["E"],
258 parsed_args["nu"],
259 parsed_args["mu_d"],
260 parsed_args["tensile_strength"],
261 parsed_args["r_min"],
262 parsed_args["r_max"],
263 parsed_args["thickness"],
264 parsed_args["rotating"],
265 parsed_args["shearvel"],
266 seed)
267 end
268
269 main()