Introduction
Introduction Statistics Contact Development Disclaimer Help
tcompress.jl - seaice-experiments - sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
LICENSE
---
tcompress.jl (8221B)
---
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 # Grain mechanical properties
102 youngs_modulus = E # Elastic modulus [Pa]
103 poissons_ratio = nu # Shear-stiffness ratio [-]
104 contact_dynamic_friction = mu_d # Coulomb-frictional coefficient [-]
105
106 # Simulation duration [s]
107 t_total = 60.0 * 60.0 * 2.0
108
109 # Temporal interval between output files
110 file_dt = 18.0
111
112 ####################################################################…
113 #### Read prior state
114 ####################################################################…
115 sim = Granular.readSimulation("rotating01-seed1/rotating01-seed1…
116 sim.id = id
117 Granular.resetTime!(sim)
118 Granular.zeroKinematics!(sim)
119 Granular.removeSimulationFiles(sim)
120
121 # Set grain mechanical properties
122 for grain in sim.grains
123 grain.youngs_modulus = youngs_modulus
124 grain.poissons_ratio = poissons_ratio
125 grain.tensile_strength = abs(tensile_strength +
126 randn()*tensile_strength_std_dev)
127 grain.shear_strength = abs(shear_strength +
128 randn()*tensile_strength_std_dev)
129 grain.contact_static_friction = contact_dynamic_friction # mu_s…
130 grain.contact_dynamic_friction = contact_dynamic_friction
131 grain.fracture_toughness = fracture_toughness
132 grain.rotating = rotating
133 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
134 end
135
136 x_max = -Inf
137 x_min = +Inf
138 y_max = -Inf
139 y_min = +Inf
140 r_max = 0.0
141 for grain in sim.grains
142 r = grain.contact_radius
143
144 if x_min > grain.lin_pos[1] - r
145 x_min = grain.lin_pos[1] - r
146 end
147
148 if x_max < grain.lin_pos[1] + r
149 x_max = grain.lin_pos[1] + r
150 end
151
152 if y_min > grain.lin_pos[2] - r
153 y_min = grain.lin_pos[2] - r
154 end
155
156 if y_max < grain.lin_pos[2] + r
157 y_max = grain.lin_pos[2] + r
158 end
159
160 if r_max < r
161 r_max = r
162 end
163 end
164
165 dx = r_max * 2.0
166 L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx]
167 n = convert(Vector{Int}, floor.(L./dx))
168 sim.ocean = Granular.createRegularOceanGrid(n, L,
169 origo=[x_min - L[1]/3.0,…
170 time=[0.], name="resized…
171 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x…
172 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y…
173
174 sim.walls[1].bc = "velocity"
175 sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total
176
177 # fix grains adjacent to walls to create no-slip BC
178 fixed_thickness = 2.0 * r_max
179 for grain in sim.grains
180 if grain.lin_pos[2] <= y_min + fixed_thickness
181 grain.fixed = true
182 grain.color = 1
183 elseif grain.lin_pos[2] >= y_max - fixed_thickness
184 grain.allow_y_acc = true
185 grain.fixed = true
186 grain.color = 1
187 end
188 end
189
190 Granular.setTimeStep!(sim)
191 Granular.setTotalTime!(sim, t_total)
192 Granular.setOutputFileInterval!(sim, file_dt)
193 render && Granular.plotGrains(sim, verbose=true)
194
195 ####################################################################…
196 #### RUN THE SIMULATION
197 ####################################################################…
198
199 # Run the consolidation experiment, and monitor top wall positio…
200 # time
201 time = Float64[]
202 compaction = Float64[]
203 effective_normal_stress = Float64[]
204 while sim.time < sim.time_total
205
206 for i=1:100
207 Granular.run!(sim, single_step=true)
208 end
209
210 append!(time, sim.time)
211 append!(compaction, sim.walls[1].pos)
212 append!(effective_normal_stress, Granular.getWallNormalS…
213
214 end
215 Granular.writeSimulation(sim)
216
217 defined_normal_stress = ones(length(effective_normal_stress)) *
218 Granular.getWallNormalStress(sim, stress_type="effective")
219 PyPlot.subplot(211)
220 PyPlot.subplots_adjust(hspace=0.0)
221 ax1 = PyPlot.gca()
222 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x l…
223 PyPlot.plot(time, compaction)
224 PyPlot.ylabel("Top wall height [m]")
225 PyPlot.subplot(212, sharex=ax1)
226 PyPlot.plot(time, defined_normal_stress)
227 PyPlot.plot(time, effective_normal_stress)
228 PyPlot.xlabel("Time [s]")
229 PyPlot.ylabel("Normal stress [Pa]")
230 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
231 PyPlot.clf()
232 DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction,
233 …
234 …
235
236 render && Granular.render(sim, trim=false, animation=true)
237 end
238
239 function main()
240 parsed_args = parse_command_line()
241 report_args(parsed_args)
242
243 seed = parsed_args["seed"]
244 run_simulation(parsed_args["id"] * "-seed" * string(seed),
245 parsed_args["E"],
246 parsed_args["nu"],
247 parsed_args["mu_d"],
248 parsed_args["tensile_strength"],
249 parsed_args["tensile_strength_std_dev"],
250 parsed_args["shear_strength"],
251 parsed_args["fracture_toughness"],
252 parsed_args["r_min"],
253 parsed_args["r_max"],
254 parsed_args["thickness"],
255 parsed_args["rotating"],
256 parsed_args["compressive_velocity"],
257 seed)
258 end
259
260 main()
261 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.