tsimulation_cons.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
---
tsimulation_cons.jl (7524B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import PyPlot
5 import ArgParse
6 using DelimitedFiles
7
8 const render = false
9
10 function parse_command_line()
11 s = ArgParse.ArgParseSettings()
12 ArgParse.@add_arg_table s begin
13 "--E"
14 help = "Young's modulus [Pa]"
15 arg_type = Float64
16 default = 2e7
17 "--nu"
18 help = "Poisson's ratio [-]"
19 arg_type = Float64
20 default = 0.285
21 "--mu_d"
22 help = "dynamic friction coefficient [-]"
23 arg_type = Float64
24 default = 0.3
25 "--gamma_n"
26 help = "contact-normal viscoity [-]"
27 arg_type = Float64
28 default = 0.0
29 "--tensile_strength"
30 help = "the maximum tensile strength [Pa]"
31 arg_type = Float64
32 default = 0.
33 "--r_min"
34 help = "minimum grain radius [m]"
35 arg_type = Float64
36 default = 5.
37 "--r_max"
38 help = "maximum grain radius [m]"
39 arg_type = Float64
40 default = 1.35e3
41 default = 50.
42 "--thickness"
43 help = "grain thickness [m]"
44 arg_type = Float64
45 default = 1.
46 "--rotating"
47 help = "allow the grains to rotate"
48 arg_type = Bool
49 default = true
50 "--shearvel"
51 help = "shear velocity [m/s]"
52 arg_type = Float64
53 default = 1.0
54 "--seed"
55 help = "seed for the pseudo RNG"
56 arg_type = Int
57 default = 1
58 "id"
59 help = "simulation id"
60 required = true
61 end
62 return ArgParse.parse_args(s)
63 end
64
65 function report_args(parsed_args)
66 println("Parsed args:")
67 for (arg,val) in parsed_args
68 println(" $arg => $val")
69 end
70 end
71
72 function run_simulation(id::String,
73 E::Float64,
74 nu::Float64,
75 mu_d::Float64,
76 gamma_n::Float64,
77 tensile_strength::Float64,
78 r_min::Float64,
79 r_max::Float64,
80 thickness::Float64,
81 rotating::Bool,
82 shearvel::Float64,
83 seed::Int)
84
85 ############################################################################
86 #### SIMULATION PARAMETERS #
87 ############################################################################
88
89 # Common simulation identifier
90 id_prefix = id
91
92 # Normal stresses for consolidation and shear [Pa]
93 N_list = [20e3]
94
95 # Simulation duration of individual steps [s]
96 t_cons = 400.0
97
98
99 ############################################################################
100 #### Step 2: Consolidate the previous output under a constant normal stress#
101 ############################################################################
102 tau_p = Float64[]
103 tau_u = Float64[]
104
105 for N in N_list
106 sim_init = Granular.createSimulation("$(id_prefix)-init")
107 sim = Granular.readSimulation(sim_init)
108
109 sim.id = "$(id_prefix)-cons-N$(N)Pa"
110 Granular.removeSimulationFiles(sim)
111
112 # Set all linear and rotational velocities to zero
113 Granular.zeroKinematics!(sim)
114
115 # Add a dynamic wall to the top which adds a normal stress downwards.
116 # The normal of this wall is downwards, and we place it at the top of
117 # the granular assemblage. Here, the fluid drag is also disabled
118 y_top = -Inf
119 for grain in sim.grains
120 grain.contact_viscosity_normal = gamma_n
121 if y_top < grain.lin_pos[2] + grain.contact_radius
122 y_top = grain.lin_pos[2] + grain.contact_radius
123 end
124 end
125 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
126 bc="normal stress",
127 normal_stress=-N)
128 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass
129 sim.walls[1].contact_viscosity_normal = 10.0*gamma_n
130 @info("Placing top wall at y=$y_top")
131
132 # Resize the grid to span the current state
133 Granular.fitGridToGrains!(sim, sim.ocean)
134
135 # Lock the grains at the very bottom so that the lower boundary is rough
136 y_bot = Inf
137 for grain in sim.grains
138 if y_bot > grain.lin_pos[2] - grain.contact_radius
139 y_bot = grain.lin_pos[2] - grain.contact_radius
140 end
141 end
142 fixed_thickness = 2.0 * r_max
143 for grain in sim.grains
144 if grain.lin_pos[2] <= fixed_thickness
145 grain.fixed = true # set x and y acceleration to zero
146 grain.color = 1
147 end
148 end
149
150 # Set current time to zero and reset output file counter
151 Granular.resetTime!(sim)
152
153 # Set the simulation time to run the consolidation for
154 Granular.setTotalTime!(sim, t_cons)
155
156 # Run the consolidation experiment, and monitor top wall position over
157 # time
158 time = Float64[]
159 compaction = Float64[]
160 effective_normal_stress = Float64[]
161 while sim.time < sim.time_total
162
163 for i=1:100
164 Granular.run!(sim, single_step=true)
165 end
166
167 append!(time, sim.time)
168 append!(compaction, sim.walls[1].pos)
169 append!(effective_normal_stress, Granular.getWallNormalStress(sim))
170
171 end
172 defined_normal_stress = ones(length(effective_normal_stress)) *
173 Granular.getWallNormalStress(sim, stress_type="effective")
174 PyPlot.subplot(211)
175 PyPlot.subplots_adjust(hspace=0.0)
176 ax1 = PyPlot.gca()
177 PyPlot.setp(ax1[:get_xticklabels](),visible=false) # Disable x labels
178 PyPlot.plot(time, compaction)
179 PyPlot.ylabel("Top wall height [m]")
180 PyPlot.subplot(212, sharex=ax1)
181 PyPlot.plot(time, defined_normal_stress)
182 PyPlot.plot(time, effective_normal_stress)
183 PyPlot.xlabel("Time [s]")
184 PyPlot.ylabel("Normal stress [Pa]")
185 PyPlot.savefig(sim.id * "-time_vs_compaction-stress.pdf")
186 PyPlot.clf()
187 writedlm(sim.id * "-data.txt", [time, compaction,
188 effective_normal_stress,
189 defined_normal_stress])
190
191 # Try to render the simulation if `pvpython` is installed on the system
192 if render
193 Granular.render(sim, trim=false)
194 end
195
196 # Save the simulation state to disk in case we need to reuse the
197 # consolidated state (e.g. different shear velocities below)
198 Granular.writeSimulation(sim)
199
200 end
201 end
202
203 function main()
204 parsed_args = parse_command_line()
205 report_args(parsed_args)
206
207 seed = parsed_args["seed"]
208 run_simulation(parsed_args["id"] * "-seed" * string(seed),
209 parsed_args["E"],
210 parsed_args["nu"],
211 parsed_args["mu_d"],
212 parsed_args["gamma_n"],
213 parsed_args["tensile_strength"],
214 parsed_args["r_min"],
215 parsed_args["r_max"],
216 parsed_args["thickness"],
217 parsed_args["rotating"],
218 parsed_args["shearvel"],
219 seed)
220 end
221
222 main()