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