Introduction
Introduction Statistics Contact Development Disclaimer Help
tinit-assemblage.jl - seaice-experiments - sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
LICENSE
---
tinit-assemblage.jl (8835B)
---
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 = 1e20
43 #default = 1285e3
44 "--r_min"
45 help = "min. grain radius [m]"
46 arg_type = Float64
47 default = 20.0
48 "--r_max"
49 help = "max. grain radius [m]"
50 arg_type = Float64
51 default = 80.0
52 "--thickness"
53 help = "grain thickness [m]"
54 arg_type = Float64
55 default = 1.
56 "--rotating"
57 help = "allow the grains to rotate"
58 arg_type = Bool
59 default = true
60 "--compressive_velocity"
61 help = "compressive velocity [m/s]"
62 arg_type = Float64
63 default = 0.1
64 "--seed"
65 help = "seed for the pseudo RNG"
66 arg_type = Int
67 default = 1
68 "id"
69 help = "simulation id"
70 required = true
71 end
72 return ArgParse.parse_args(s)
73 end
74
75 function report_args(parsed_args)
76 println("Parsed args:")
77 for (arg,val) in parsed_args
78 println(" $arg => $val")
79 end
80 end
81
82 function run_simulation(id::String,
83 E::Float64,
84 nu::Float64,
85 mu_d::Float64,
86 tensile_strength::Float64,
87 tensile_strength_std_dev::Float64,
88 shear_strength::Float64,
89 fracture_toughness::Float64,
90 r_min::Float64,
91 r_max::Float64,
92 thickness::Float64,
93 rotating::Bool,
94 compressive_velocity::Float64,
95 seed::Int)
96
97 ####################################################################…
98 #### SIMULATION PARAMETERS
99 ####################################################################…
100
101 # Common simulation identifier
102 id_prefix = id
103
104 # Grain mechanical properties
105 youngs_modulus = E # Elastic modulus [Pa]
106 poissons_ratio = nu # Shear-stiffness ratio [-]
107 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
108
109 # Simulation duration [s]
110 t_total = 60.0 * 60.0 * 4.0
111
112 # Temporal interval between output files
113 file_dt = 18.0
114
115 # Misc parameters
116 water_density = 1000.
117
118 # World size [m]
119 L = [5e3, 2e4, 10.0]
120
121 Random.seed!(seed)
122
123 ####################################################################…
124 #### Create ice floes
125 ####################################################################…
126 sim = Granular.createSimulation("$(id_prefix)")
127 Granular.removeSimulationFiles(sim)
128
129 # Create box edges of ice floes with r_min radius, moving inward
130 r_walls = r_min
131 #= for x in LinRange(0.0, L[1] - r_walls, Int(ceil(L[1]/(r_walls*2))…
132 #= Granular.addGrainCylindrical!(sim, [x, r_min], =#
133 #= r_walls, thickness, fixed=true, …
134 #= lin_vel=[0.0, compressive_veloci…
135 #= verbose=false) =#
136 #= Granular.addGrainCylindrical!(sim, [x, L[2] - r_min], =#
137 #= r_walls, thickness, fixed=true, …
138 #= lin_vel=[0.0, -compressive_veloc…
139 #= verbose=false) =#
140 #= end =#
141 #for y in LinRange(r_walls, L[2] - r_walls, Int(ceil(L[2]/(r_walls*2…
142 # Granular.addGrainCylindrical!(sim, [r_min, y],
143 # r_walls, thickness, fixed=true,
144 # #lin_vel=[compressive_velocity/2.,…
145 # verbose=false)
146 # Granular.addGrainCylindrical!(sim, [L[1] - r_min, y],
147 # r_walls, thickness, fixed=true,
148 # lin_vel=[-compressive_velocity/2.,…
149 # verbose=false)
150 #end
151
152 # Create a grid for contact searching spanning the extent of the gra…
153 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2]
154 sim.ocean = Granular.createRegularOceanGrid(n, L)
155 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x…
156 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y…
157
158 # Add unconstrained ice floes
159 #= Granular.regularPacking!(sim, =#
160 #= n, =#
161 #= r_min, r_max, =#
162 #= seed=seed) =#
163 Granular.irregularPacking!(sim,
164 radius_max=r_max, radius_min=r_min,
165 thickness=thickness,
166 #binary_radius_search=true,
167 seed=seed)
168
169 # Set grain mechanical properties
170 for grain in sim.grains
171 grain.youngs_modulus = youngs_modulus
172 grain.poissons_ratio = poissons_ratio
173 grain.tensile_strength = abs(tensile_strength +
174 randn()*tensile_strength_std_dev)
175 grain.shear_strength = abs(shear_strength +
176 randn()*tensile_strength_std_dev)
177 grain.contact_static_friction = contact_dynamic_friction # mu_s…
178 grain.contact_dynamic_friction = contact_dynamic_friction
179 grain.fracture_toughness = fracture_toughness
180 grain.rotating = rotating
181 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
182 end
183
184 # Create GSD plot to signify that simulation is complete
185 Granular.writeSimulation(sim)
186 render && Granular.plotGrainSizeDistribution(sim)
187
188 # Add a dynamic wall to the top which adds a normal stress downw…
189 # The normal of this wall is downwards, and we place it at the t…
190 # the granular assemblage. Here, the fluid drag is also disabled
191 y_top = -Inf
192 for grain in sim.grains
193 if y_top < grain.lin_pos[2] + grain.contact_radius
194 y_top = grain.lin_pos[2] + grain.contact_radius
195 end
196 end
197 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
198 …
199 …
200 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain …
201 @info("Placing top wall at y=$y_top")
202
203 # Resize the grid to span the current state
204 Granular.fitGridToGrains!(sim, sim.ocean)
205
206 Granular.setTimeStep!(sim)
207 Granular.setTotalTime!(sim, t_total)
208 Granular.setOutputFileInterval!(sim, file_dt)
209 render && Granular.plotGrains(sim, verbose=true)
210
211 ####################################################################…
212 #### RUN THE SIMULATION
213 ####################################################################…
214
215 Granular.run!(sim)
216 Granular.writeSimulation(sim)
217
218 render && Granular.render(sim, trim=false, animation=true)
219 end
220
221 function main()
222 parsed_args = parse_command_line()
223 report_args(parsed_args)
224
225 seed = parsed_args["seed"]
226 run_simulation(parsed_args["id"] * "-seed" * string(seed),
227 parsed_args["E"],
228 parsed_args["nu"],
229 parsed_args["mu_d"],
230 parsed_args["tensile_strength"],
231 parsed_args["tensile_strength_std_dev"],
232 parsed_args["shear_strength"],
233 parsed_args["fracture_toughness"],
234 parsed_args["r_min"],
235 parsed_args["r_max"],
236 parsed_args["thickness"],
237 parsed_args["rotating"],
238 parsed_args["compressive_velocity"],
239 seed)
240 end
241
242 main()
243 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.