tcompress.jl - seaice-exp-faulting - simulate faulting angles in sea ice under uniaxial compression
(HTM) git clone git://src.adamsgaard.dk/seaice-exp-faulting
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tcompress.jl (8238B)
---
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 = 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.681.jld2")
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 = mu_d
130 grain.contact_dynamic_friction = contact_dynamic_friction
131 grain.fracture_toughness = fracture_toughness
132 grain.rotating = rotating
133 grain.ocean_drag_coeff_horiz = 1.6e-4
134 grain.ocean_drag_coeff_vert = 0.14
135 end
136
137 x_max = -Inf
138 x_min = +Inf
139 y_max = -Inf
140 y_min = +Inf
141 r_max = 0.0
142 for grain in sim.grains
143 r = grain.contact_radius
144
145 if x_min > grain.lin_pos[1] - r
146 x_min = grain.lin_pos[1] - r
147 end
148
149 if x_max < grain.lin_pos[1] + r
150 x_max = grain.lin_pos[1] + r
151 end
152
153 if y_min > grain.lin_pos[2] - r
154 y_min = grain.lin_pos[2] - r
155 end
156
157 if y_max < grain.lin_pos[2] + r
158 y_max = grain.lin_pos[2] + r
159 end
160
161 if r_max < r
162 r_max = r
163 end
164 end
165
166 dx = r_max * 2.0
167 L::Vector{Float64} = [(x_max - x_min)*3.0, y_max - y_min, dx]
168 n = convert(Vector{Int}, floor.(L./dx))
169 sim.ocean = Granular.createRegularOceanGrid(n, L,
170 origo=[x_min - L[1]/3.0, y_min],
171 time=[0.], name="resized")
172 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
173 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
174
175 sim.walls[1].bc = "velocity"
176 sim.walls[1].vel = -(y_max - y_min) * 0.2 / t_total
177
178 # fix grains adjacent to walls to create no-slip BC
179 fixed_thickness = 2.0 * r_max
180 for grain in sim.grains
181 if grain.lin_pos[2] <= y_min + fixed_thickness
182 grain.fixed = true
183 grain.color = 1
184 elseif grain.lin_pos[2] >= y_max - fixed_thickness
185 grain.allow_y_acc = true
186 grain.fixed = true
187 grain.color = 1
188 end
189 end
190
191 Granular.setTimeStep!(sim)
192 Granular.setTotalTime!(sim, t_total)
193 Granular.setOutputFileInterval!(sim, file_dt)
194 render && Granular.plotGrains(sim, verbose=true)
195
196 ############################################################################
197 #### RUN THE SIMULATION
198 ############################################################################
199
200 # Run the consolidation experiment, and monitor top wall position over
201 # time
202 time = Float64[]
203 compaction = Float64[]
204 effective_normal_stress = Float64[]
205 while sim.time < sim.time_total
206
207 for i=1:100
208 Granular.run!(sim, single_step=true)
209 end
210
211 append!(time, sim.time)
212 append!(compaction, sim.walls[1].pos)
213 append!(effective_normal_stress, Granular.getWallNormalStress(sim))
214
215 end
216 Granular.writeSimulation(sim)
217
218 defined_normal_stress = ones(length(effective_normal_stress)) *
219 Granular.getWallNormalStress(sim, stress_type="effective")
220 PyPlot.subplot(211)
221 PyPlot.subplots_adjust(hspace=0.0)
222 ax1 = PyPlot.gca()
223 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
224 PyPlot.plot(time, compaction)
225 PyPlot.ylabel("Top wall height [m]")
226 PyPlot.subplot(212, sharex=ax1)
227 PyPlot.plot(time, defined_normal_stress)
228 PyPlot.plot(time, effective_normal_stress)
229 PyPlot.xlabel("Time [s]")
230 PyPlot.ylabel("Normal stress [Pa]")
231 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
232 PyPlot.clf()
233 DelimitedFiles.writedlm(sim.id * "-data.txt", [time, compaction,
234 effective_normal_stress,
235 defined_normal_stress])
236
237 render && Granular.render(sim, trim=false, animation=true)
238 end
239
240 function main()
241 parsed_args = parse_command_line()
242 report_args(parsed_args)
243
244 seed = parsed_args["seed"]
245 run_simulation(parsed_args["id"] * "-seed" * string(seed),
246 parsed_args["E"],
247 parsed_args["nu"],
248 parsed_args["mu_d"],
249 parsed_args["tensile_strength"],
250 parsed_args["tensile_strength_std_dev"],
251 parsed_args["shear_strength"],
252 parsed_args["fracture_toughness"],
253 parsed_args["r_min"],
254 parsed_args["r_max"],
255 parsed_args["thickness"],
256 parsed_args["rotating"],
257 parsed_args["compressive_velocity"],
258 seed)
259 end
260
261 main()
262 end