tridging_bulk_simulation.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
---
tridging_bulk_simulation.jl (9545B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import PyPlot
5 import ArgParse
6 using Random
7 import DelimitedFiles
8
9 let
10 render = false
11
12 function parse_command_line()
13 s = ArgParse.ArgParseSettings()
14 ArgParse.@add_arg_table s begin
15 "--E"
16 help = "Young's modulus [Pa]"
17 arg_type = Float64
18 default = 2e7
19 "--nu"
20 help = "Poisson's ratio [-]"
21 arg_type = Float64
22 default = 0.285
23 "--mu_d"
24 help = "dynamic friction coefficient [-]"
25 arg_type = Float64
26 default = 0.3
27 "--tensile_strength"
28 help = "the maximum tensile strength [Pa]"
29 arg_type = Float64
30 default = 400e3
31 "--tensile_strength_std_dev"
32 help = "standard deviation of the maximum tensile strength [Pa]"
33 arg_type = Float64
34 default = 0.0
35 "--shear_strength"
36 help = "the maximum shear strength [Pa]"
37 arg_type = Float64
38 default = 200e3
39 "--fracture_toughness"
40 help = "fracture toughness [Pa*m^{1/2}]"
41 arg_type = Float64
42 default = 1285e3
43 "--r_min"
44 help = "min. grain radius [m]"
45 arg_type = Float64
46 default = 20.0
47 "--r_max"
48 help = "max. grain radius [m]"
49 arg_type = Float64
50 default = 200.0
51 "--thickness"
52 help = "grain thickness [m]"
53 arg_type = Float64
54 default = 1.
55 "--rotating"
56 help = "allow the grains to rotate"
57 arg_type = Bool
58 default = true
59 "--compressive_velocity"
60 help = "compressive velocity [m/s]"
61 arg_type = Float64
62 default = 0.1
63 "--seed"
64 help = "seed for the pseudo RNG"
65 arg_type = Int
66 default = 1
67 "id"
68 help = "simulation id"
69 required = true
70 end
71 return ArgParse.parse_args(s)
72 end
73
74 function report_args(parsed_args)
75 println("Parsed args:")
76 for (arg,val) in parsed_args
77 println(" $arg => $val")
78 end
79 end
80
81 function run_simulation(id::String,
82 E::Float64,
83 nu::Float64,
84 mu_d::Float64,
85 tensile_strength::Float64,
86 tensile_strength_std_dev::Float64,
87 shear_strength::Float64,
88 fracture_toughness::Float64,
89 r_min::Float64,
90 r_max::Float64,
91 thickness::Float64,
92 rotating::Bool,
93 compressive_velocity::Float64,
94 seed::Int)
95
96 ############################################################################
97 #### SIMULATION PARAMETERS
98 ############################################################################
99
100 # Common simulation identifier
101 id_prefix = id
102
103 # Grain mechanical properties
104 youngs_modulus = E # Elastic modulus [Pa]
105 poissons_ratio = nu # Shear-stiffness ratio [-]
106 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
107
108 # Simulation duration [s]
109 t_total = 60.0 * 60.0
110
111 # Temporal interval between output files
112 file_dt = t_total/200
113
114 # Misc parameters
115 water_density = 1000.
116
117 # World size [m]
118 L = [1e3, 1e3, 10.0]
119
120 Random.seed!(seed)
121
122 ############################################################################
123 #### Create ice floes
124 ############################################################################
125 sim = Granular.createSimulation("$(id_prefix)")
126 Granular.removeSimulationFiles(sim)
127
128 # Create box edges of ice floes with r_min radius, moving inward
129 r_walls = r_min
130 #= for x in LinRange(0.0, L[1] - r_walls, Int(ceil(L[1]/(r_walls*2)))) =#
131 #= Granular.addGrainCylindrical!(sim, [x, r_min], =#
132 #= r_walls, thickness, fixed=true, =#
133 #= lin_vel=[0.0, compressive_velocity/2.], =#
134 #= verbose=false) =#
135 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =#
136 #= r_walls, thickness, fixed=true, =#
137 #= lin_vel=[0.0, -compressive_velocity/2.], =#
138 #= verbose=false) =#
139 #= end =#
140 for y in LinRange(r_walls, L[2] - r_walls, Int(ceil(L[2]/(r_walls*2))))
141 Granular.addGrainCylindrical!(sim, [r_min, y],
142 r_walls, thickness, fixed=true,
143 #lin_vel=[compressive_velocity/2., 0.0],
144 verbose=false)
145 Granular.addGrainCylindrical!(sim, [L[1] - r_min, y],
146 r_walls, thickness, fixed=true,
147 lin_vel=[-compressive_velocity/2., 0.0],
148 verbose=false)
149 end
150
151 # Create a grid for contact searching spanning the extent of the grains
152 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2]
153 sim.ocean = Granular.createRegularOceanGrid(n, L)
154 Granular.setGridBoundaryConditions!(sim.ocean, "inactive", "-x +x")
155 Granular.setGridBoundaryConditions!(sim.ocean, "periodic", "-y +y")
156
157 # Add unconstrained ice floes
158 #= Granular.regularPacking!(sim, =#
159 #= n, =#
160 #= r_min, r_max, =#
161 #= seed=seed) =#
162 Granular.irregularPacking!(sim,
163 radius_max=r_max, radius_min=r_min,
164 thickness=thickness,
165 #binary_radius_search=true,
166 seed=seed)
167
168 # Set grain mechanical properties
169 for grain in sim.grains
170 grain.youngs_modulus = youngs_modulus
171 grain.poissons_ratio = poissons_ratio
172 grain.tensile_strength = abs(tensile_strength +
173 randn()*tensile_strength_std_dev)
174 grain.shear_strength = abs(shear_strength +
175 randn()*tensile_strength_std_dev)
176 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
177 grain.contact_dynamic_friction = contact_dynamic_friction
178 grain.fracture_toughness = fracture_toughness
179 grain.rotating = rotating
180 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
181 end
182
183 # Create GSD plot to signify that simulation is complete
184 Granular.writeSimulation(sim)
185 render && Granular.plotGrainSizeDistribution(sim)
186
187 Granular.setTimeStep!(sim)
188 Granular.setTotalTime!(sim, t_total)
189 Granular.setOutputFileInterval!(sim, file_dt)
190 render && Granular.plotGrains(sim, verbose=true)
191
192
193 ############################################################################
194 #### RUN THE SIMULATION
195 ############################################################################
196 #= theta = 0.0; submerged_volume = 0.0 =#
197 time = Float64[]
198 N = Float64[] # normal stress on leftmost fixed particles
199 τ = Float64[] # shear stress on leftmost fixed particles
200 A = sim.grains[1].thickness * L[2]
201 while sim.time < sim.time_total
202
203 append!(time, sim.time)
204
205 # Record cumulative force on fixed particles at the +x boundary
206 normal_force = 0.0
207 tangential_force = 0.0
208 for grain in sim.grains
209 if grain.fixed && grain.lin_pos[1] > 0.2*L[1]
210 normal_force += grain.force[1]
211 tangential_force += grain.force[2]
212 end
213 end
214 append!(N, normal_force/A) # compressive = positive
215 append!(τ, tangential_force/A)
216
217 # Run for 100 time steps before measuring bulk forces
218 for i=1:100
219 Granular.run!(sim, single_step=true)
220 end
221
222 end
223
224 # Plot normal stress and shear stress vs time
225 DelimitedFiles.writedlm(id_prefix * "-data.txt", [time, N, τ])
226 PyPlot.subplot(211)
227 PyPlot.subplots_adjust(hspace=0.0)
228 ax1 = PyPlot.gca()
229 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
230 PyPlot.plot(time, N/1e3)
231 PyPlot.ylabel("Compressive stress [kPa]")
232 PyPlot.subplot(212, sharex=ax1)
233 PyPlot.plot(time, τ/1e3)
234 PyPlot.xlabel("Time [s]")
235 PyPlot.ylabel("Shear stress [kPa]")
236 PyPlot.tight_layout()
237 PyPlot.savefig(sim.id * "-time_vs_stress.pdf")
238 PyPlot.clf()
239
240 render && Granular.render(sim, trim=false, animation=true)
241 end
242
243 function main()
244 parsed_args = parse_command_line()
245 report_args(parsed_args)
246
247 seed = parsed_args["seed"]
248 run_simulation(parsed_args["id"] * "-seed" * string(seed),
249 parsed_args["E"],
250 parsed_args["nu"],
251 parsed_args["mu_d"],
252 parsed_args["tensile_strength"],
253 parsed_args["tensile_strength_std_dev"],
254 parsed_args["shear_strength"],
255 parsed_args["fracture_toughness"],
256 parsed_args["r_min"],
257 parsed_args["r_max"],
258 parsed_args["thickness"],
259 parsed_args["rotating"],
260 parsed_args["compressive_velocity"],
261 seed)
262 end
263
264 main()
265 end