Introduction
Introduction Statistics Contact Development Disclaimer Help
tridging_bulk_simulation.jl - seaice-experiments - sea ice experiments using Gr…
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
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 […
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_veloci…
134 #= verbose=false) =#
135 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =#
136 #= r_walls, thickness, fixed=true, …
137 #= lin_vel=[0.0, -compressive_veloc…
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., …
144 verbose=false)
145 Granular.addGrainCylindrical!(sim, [L[1] - r_min, y],
146 r_walls, thickness, fixed=true,
147 lin_vel=[-compressive_velocity/2., …
148 verbose=false)
149 end
150
151 # Create a grid for contact searching spanning the extent of the gra…
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…
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
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.