tinit-assemblage.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
---
tinit-assemblage.jl (7396B)
---
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 # 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
132 # Create a grid for contact searching spanning the extent of the grains
133 n = [Int(ceil(L[1]/r_max/2)), Int(ceil(L[2]/r_max/2)), 2]
134 sim.ocean = Granular.createRegularOceanGrid(n, L)
135 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-x +x")
136 Granular.setGridBoundaryConditions!(sim.ocean, "impermeable", "-y +y")
137
138 # Add unconstrained ice floes
139 Granular.irregularPacking!(sim,
140 radius_max=r_max, radius_min=r_min,
141 thickness=thickness,
142 #binary_radius_search=true,
143 seed=seed)
144
145 # Set grain mechanical properties
146 for grain in sim.grains
147 grain.youngs_modulus = youngs_modulus
148 grain.poissons_ratio = poissons_ratio
149 grain.tensile_strength = abs(tensile_strength +
150 randn()*tensile_strength_std_dev)
151 grain.shear_strength = abs(shear_strength +
152 randn()*tensile_strength_std_dev)
153 grain.contact_static_friction = contact_dynamic_friction # mu_s = mu_d
154 grain.contact_dynamic_friction = contact_dynamic_friction
155 grain.fracture_toughness = fracture_toughness
156 grain.rotating = rotating
157 grain.ocean_drag_coeff_horiz = 0.0 # don't include drag on ends
158 end
159
160 # Create GSD plot to signify that simulation is complete
161 Granular.writeSimulation(sim)
162 render && Granular.plotGrainSizeDistribution(sim)
163
164 # Add a dynamic wall to the top which adds a normal stress downwards.
165 # The normal of this wall is downwards, and we place it at the top of
166 # the granular assemblage. Here, the fluid drag is also disabled
167 y_top = -Inf
168 for grain in sim.grains
169 if y_top < grain.lin_pos[2] + grain.contact_radius
170 y_top = grain.lin_pos[2] + grain.contact_radius
171 end
172 end
173 Granular.addWallLinearFrictionless!(sim, [0., 1.], y_top,
174 bc="normal stress",
175 normal_stress=-20e3)
176 sim.walls[1].mass *= 0.1 # set wall mass to 10% of total grain mass
177 @info("Placing top wall at y=$y_top")
178
179 # Resize the grid to span the current state
180 Granular.fitGridToGrains!(sim, sim.ocean)
181
182 Granular.setTimeStep!(sim)
183 Granular.setTotalTime!(sim, t_total)
184 Granular.setOutputFileInterval!(sim, file_dt)
185 render && Granular.plotGrains(sim, verbose=true)
186
187 ############################################################################
188 #### RUN THE SIMULATION
189 ############################################################################
190
191 Granular.run!(sim)
192 Granular.writeSimulation(sim)
193
194 render && Granular.render(sim, trim=false, animation=true)
195 end
196
197 function main()
198 parsed_args = parse_command_line()
199 report_args(parsed_args)
200
201 seed = parsed_args["seed"]
202 run_simulation(parsed_args["id"] * "-seed" * string(seed),
203 parsed_args["E"],
204 parsed_args["nu"],
205 parsed_args["mu_d"],
206 parsed_args["tensile_strength"],
207 parsed_args["tensile_strength_std_dev"],
208 parsed_args["shear_strength"],
209 parsed_args["fracture_toughness"],
210 parsed_args["r_min"],
211 parsed_args["r_max"],
212 parsed_args["thickness"],
213 parsed_args["rotating"],
214 parsed_args["compressive_velocity"],
215 seed)
216 end
217
218 main()
219 end