Introduction
Introduction Statistics Contact Development Disclaimer Help
tsimulation_cons.jl - seaice-experiments - sea ice experiments using Granular.jl
git clone git://src.adamsgaard.dk/seaice-experiments
Log
Files
Refs
README
LICENSE
---
tsimulation_cons.jl (7577B)
---
1 #/usr/bin/env julia
2 ENV["MPLBACKEND"] = "Agg"
3 import Granular
4 import JLD
5 import PyPlot
6 import ArgParse
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 const id_prefix = id
91
92 # Normal stresses for consolidation and shear [Pa]
93 const N_list = [5e3, 10e3, 15e3, 20e3, 30e3, 40e3, 80e3, 160e3]
94
95 # Simulation duration of individual steps [s]
96 const t_cons = 400.0
97
98
99 ####################################################################…
100 #### Step 2: Consolidate the previous output under a constant normal…
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 downw…
116 # The normal of this wall is downwards, and we place it at the t…
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 …
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 …
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 const 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 positio…
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.getWallNormalStres…
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 l…
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…
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()
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.