 |
tgrain.jl - Granular.jl - Julia package for granular dynamics simulation |
 |
git clone git://src.adamsgaard.dk/Granular.jl |
 |
Log |
|
Files |
|
Refs |
|
README |
|
LICENSE |
|
--- |
|
tgrain.jl (39136B) |
|
--- |
|
1 ## Manage grains in the model |
|
2 |
|
3 using Test |
|
4 using LinearAlgebra |
|
5 |
|
6 export addGrainCylindrical! |
|
7 """ |
|
8 function addGrainCylindrical!(simulation, lin_pos, contact_radius, |
|
9 thickness[, areal_radius, lin_vel, l… |
|
10 force, ang_pos, ang_vel, ang_acc, to… |
|
11 density, contact_stiffness_normal, |
|
12 contact_stiffness_tangential, |
|
13 contact_viscosity_normal, |
|
14 contact_viscosity_tangential, |
|
15 contact_static_friction, |
|
16 contact_dynamic_friction, |
|
17 youngs_modulus, poissons_ratio, |
|
18 tensile_strength, shear_strength, |
|
19 strength_heal_rate, |
|
20 fracture_toughness, |
|
21 ocean_drag_coeff_vert, |
|
22 ocean_drag_coeff_horiz, |
|
23 atmosphere_drag_coeff_vert, |
|
24 atmosphere_drag_coeff_horiz, |
|
25 pressure, fixed, |
|
26 allow_x_acc, allow_y_acc, allow_z_ac… |
|
27 rotating, enabled, verbose, |
|
28 ocean_grid_pos, atmosphere_grid_pos, |
|
29 n_contact, granular_stress, ocean_st… |
|
30 atmosphere_stress, |
|
31 thermal_energy, |
|
32 color]) |
|
33 |
|
34 Creates and adds a cylindrical grain to a simulation. Most of the argume… |
|
35 are optional, and come with default values. The only required arguments… |
|
36 `simulation`, `lin_pos`, `contact_radius`, and `thickness`. |
|
37 |
|
38 # Arguments |
|
39 * `simulation::Simulation`: the simulation object where the grain should… |
|
40 added to. |
|
41 * `lin_pos::Vector{Float64}`: linear position of grain center [m]. If a |
|
42 two-component vector is used, the values will be mapped to *x* and *… |
|
43 the *z* component will be set to zero. |
|
44 * `contact_radius::Float64`: grain radius for granular interaction [m]. |
|
45 * `thickness::Float64`: grain thickness [m]. |
|
46 * `areal_radius = false`: grain radius for determining sea-ice concentra… |
|
47 [m]. |
|
48 * `lin_vel::Vector{Float64} = [0., 0., 0.]`: linear velocity [m/s]. If a |
|
49 two-component vector is used, the values will be mapped to *x* and *… |
|
50 the *z* component will be set to zero. |
|
51 * `lin_acc::Vector{Float64} = [0., 0., 0.]`: linear acceleration [m/s^2]… |
|
52 two-component vector is used, the values will be mapped to *x* and *… |
|
53 the *z* component will be set to zero. |
|
54 * `force::Vector{Float64} = [0., 0., 0.]`: linear force balance [N]. If a |
|
55 two-component vector is used, the values will be mapped to *x* and *… |
|
56 the *z* component will be set to zero. |
|
57 * `ang_pos::Float64 = [0., 0., 0.]`: angular position around its center … |
|
58 axis [rad]. If a scalar is used, the value will be mapped to *z*, an… |
|
59 *x* and *y* components will be set to zero. |
|
60 * `ang_vel::Float64 = [0., 0., 0.]`: angular velocity around its center … |
|
61 axis [rad/s]. If a scalar is used, the value will be mapped to *z*, … |
|
62 *x* and *y* components will be set to zero. |
|
63 * `ang_acc::Float64 = [0., 0., 0.]`: angular acceleration around its cen… |
|
64 vertical axis [rad/s^2]. If a scalar is used, the value will be mapp… |
|
65 *z*, and the *x* and *y* components will be set to zero. |
|
66 * `torque::Float64 = [0., 0., 0.]`: torque around its center vertical ax… |
|
67 [N*m]. If a scalar is used, the value will be mapped to *z*, and the… |
|
68 *y* components will be set to zero. |
|
69 * `density::Float64 = 934.`: grain mean density [kg/m^3]. |
|
70 * `contact_stiffness_normal::Float64 = 1e7`: contact-normal stiffness [N… |
|
71 overridden if `youngs_modulus` is set to a positive value. |
|
72 * `contact_stiffness_tangential::Float64 = 0.`: contact-tangential stiff… |
|
73 [N/m]; overridden if `youngs_modulus` is set to a positive value. |
|
74 * `contact_viscosity_normal::Float64 = 0.`: contact-normal viscosity [N/… |
|
75 * `contact_viscosity_tangential::Float64 = 0.`: contact-tangential visco… |
|
76 [N/m/s]. |
|
77 * `contact_static_friction::Float64 = 0.4`: contact static Coulomb frict… |
|
78 coefficient [-]. |
|
79 * `contact_dynamic_friction::Float64 = 0.4`: contact dynamic Coulomb fri… |
|
80 coefficient [-]. |
|
81 * `youngs_modulus::Float64 = 2e7`: elastic modulus [Pa]; overrides any v… |
|
82 set for `contact_stiffness_normal`. |
|
83 * `poissons_ratio::Float64 = 0.185`: Poisson's ratio, used to determine … |
|
84 contact-tangential stiffness from `youngs_modulus` [-]. |
|
85 * `tensile_strength::Float64 = 0.`: contact-tensile (cohesive) bond stre… |
|
86 [Pa]. |
|
87 * `shear_strength::Float64 = 0.`: shear strength of bonded contacts [Pa]. |
|
88 * `strength_heal_rate::Float64 = 0.`: rate at which contact bond |
|
89 strength is obtained [Pa/s]. |
|
90 * `fracture_toughness::Float64 = 0.`: fracture toughness which influence… |
|
91 maximum compressive strength on granular contact [m^{1/2}*Pa]. A val… |
|
92 of 1.285e6 m^{1/2}*Pa is used for sea ice by Hopkins 2004. |
|
93 * `ocean_drag_coeff_vert::Float64 = 0.85`: vertical drag coefficient for… |
|
94 against grain sides [-]. |
|
95 * `ocean_drag_coeff_horiz::Float64 = 5e-4`: horizontal drag coefficient … |
|
96 ocean against grain bottom [-]. |
|
97 * `atmosphere_drag_coeff_vert::Float64 = 0.4`: vertical drag coefficient… |
|
98 atmosphere against grain sides [-]. |
|
99 * `atmosphere_drag_coeff_horiz::Float64 = 2.5e-4`: horizontal drag coeff… |
|
100 for atmosphere against grain bottom [-]. |
|
101 * `pressure::Float64 = 0.`: current compressive stress on grain [Pa]. |
|
102 * `fixed::Bool = false`: grain is fixed to a constant velocity (e.g. zer… |
|
103 * `allow_x_acc::Bool = false`: override `fixed` along `x`. |
|
104 * `allow_y_acc::Bool = false`: override `fixed` along `y`. |
|
105 * `allow_z_acc::Bool = false`: override `fixed` along `z`. |
|
106 * `rotating::Bool = true`: grain is allowed to rotate. |
|
107 * `enabled::Bool = true`: grain interacts with other grains. |
|
108 * `verbose::Bool = true`: display diagnostic information during the func… |
|
109 call. |
|
110 * `ocean_grid_pos::Array{Int, 1} = [0, 0]`: position of grain in the oce… |
|
111 grid. |
|
112 * `atmosphere_grid_pos::Array{Int, 1} = [0, 0]`: position of grain in the |
|
113 atmosphere grid. |
|
114 * `n_contacts::Int = 0`: number of contacts with other grains. |
|
115 * `granular_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on… |
|
116 from granular interactions [Pa]. |
|
117 * `ocean_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on gr… |
|
118 ocean drag [Pa]. |
|
119 * `atmosphere_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress … |
|
120 from atmosphere drag [Pa]. |
|
121 * `thermal_energy::Float64 = 0.0`: thermal energy of grain [J]. |
|
122 * `color::Int=0`: type number, usually used for associating a color to t… |
|
123 during visualization. |
|
124 |
|
125 # Examples |
|
126 The most basic example adds a new grain to the simulation `sim`, with a |
|
127 center at `[1., 2., 0.]`, a radius of `1.` meter, and a thickness of `0.… |
|
128 meter: |
|
129 |
|
130 ```julia |
|
131 Granular.addGrainCylindrical!(sim, [1., 2.], 1., .5) |
|
132 ``` |
|
133 Note that the *z* component is set to zero if a two-component vector is … |
|
134 |
|
135 The following example will create a grain with tensile and shear strengt… |
|
136 velocity of 0.5 m/s towards -x: |
|
137 |
|
138 ```julia |
|
139 Granular.addGrainCylindrical!(sim, [4., 2.], 1., .5, |
|
140 tensile_strength = 200e3, |
|
141 shear_strength = 100e3, |
|
142 lin_vel = [-.5, 0.]) |
|
143 ``` |
|
144 |
|
145 Fixed grains are useful for creating walls or coasts, and loops are usef… |
|
146 for creating regular arrangements: |
|
147 |
|
148 ```julia |
|
149 for i=1:5 |
|
150 Granular.addGrainCylindrical!(sim, [i*2., 0., 3.], 1., .5, fixed=tru… |
|
151 end |
|
152 ``` |
|
153 """ |
|
154 function addGrainCylindrical!(simulation::Simulation, |
|
155 lin_pos::Vector{Float64}, |
|
156 contact_radius::Float64, |
|
157 thickness::Float64; |
|
158 areal_radius = false, |
|
159 lin_vel::Vector{Float64} = [0., 0., 0.], |
|
160 lin_acc::Vector{Float64} = [0., 0., 0.], |
|
161 force::Vector{Float64} = [0., 0., 0.], |
|
162 ang_pos::Vector{Float64} = [0., 0., 0.], |
|
163 ang_vel::Vector{Float64} = [0., 0., 0.], |
|
164 ang_acc::Vector{Float64} = [0., 0., 0.], |
|
165 torque::Vector{Float64} = [0., 0., 0.], |
|
166 density::Float64 = 934., |
|
167 contact_stiffness_normal::Float64 = 1e7, |
|
168 contact_stiffness_tangential::Float64 = … |
|
169 contact_viscosity_normal::Float64 = 0., |
|
170 contact_viscosity_tangential::Float64 = … |
|
171 contact_static_friction::Float64 = 0.4, |
|
172 contact_dynamic_friction::Float64 = 0.4, |
|
173 youngs_modulus::Float64 = 2e7, |
|
174 poissons_ratio::Float64 = 0.185, # Hopk… |
|
175 tensile_strength::Float64 = 0., |
|
176 shear_strength::Float64 = 0., |
|
177 strength_heal_rate::Float64 = Inf, |
|
178 fracture_toughness::Float64 = 0., |
|
179 ocean_drag_coeff_vert::Float64 = 0.85, #… |
|
180 ocean_drag_coeff_horiz::Float64 = 5e-4, … |
|
181 atmosphere_drag_coeff_vert::Float64 = 0.… |
|
182 atmosphere_drag_coeff_horiz::Float64 = 2… |
|
183 pressure::Float64 = 0., |
|
184 fixed::Bool = false, |
|
185 allow_x_acc::Bool = false, |
|
186 allow_y_acc::Bool = false, |
|
187 allow_z_acc::Bool = false, |
|
188 rotating::Bool = true, |
|
189 enabled::Bool = true, |
|
190 verbose::Bool = true, |
|
191 ocean_grid_pos::Array{Int, 1} = [0, 0], |
|
192 atmosphere_grid_pos::Array{Int, 1} = [0,… |
|
193 n_contacts::Int = 0, |
|
194 granular_stress::Vector{Float64} = [0., … |
|
195 ocean_stress::Vector{Float64} = [0., 0.,… |
|
196 atmosphere_stress::Vector{Float64} = [0.… |
|
197 thermal_energy::Float64 = 0., |
|
198 color::Int = 0) |
|
199 |
|
200 # Check input values |
|
201 if length(lin_pos) != 3 |
|
202 lin_pos = vecTo3d(lin_pos) |
|
203 end |
|
204 if length(lin_vel) != 3 |
|
205 lin_vel = vecTo3d(lin_vel) |
|
206 end |
|
207 if length(lin_acc) != 3 |
|
208 lin_acc = vecTo3d(lin_acc) |
|
209 end |
|
210 if length(force) != 3 |
|
211 force = vecTo3d(force) |
|
212 end |
|
213 if length(ang_pos) != 3 |
|
214 ang_pos = vecTo3d(ang_pos) |
|
215 end |
|
216 if length(ang_vel) != 3 |
|
217 ang_vel = vecTo3d(ang_vel) |
|
218 end |
|
219 if length(ang_acc) != 3 |
|
220 ang_acc = vecTo3d(ang_acc) |
|
221 end |
|
222 if length(torque) != 3 |
|
223 torque = vecTo3d(torque) |
|
224 end |
|
225 if length(granular_stress) != 3 |
|
226 granular_stress = vecTo3d(granular_stress) |
|
227 end |
|
228 if length(ocean_stress) != 3 |
|
229 ocean_stress = vecTo3d(ocean_stress) |
|
230 end |
|
231 if length(atmosphere_stress) != 3 |
|
232 atmosphere_stress = vecTo3d(atmosphere_stress) |
|
233 end |
|
234 if contact_radius <= 0.0 |
|
235 error("Radius must be greater than 0.0 " * |
|
236 "(radius = $contact_radius m)") |
|
237 end |
|
238 if density <= 0.0 |
|
239 error("Density must be greater than 0.0 " * |
|
240 "(density = $density kg/m^3)") |
|
241 end |
|
242 |
|
243 if !areal_radius |
|
244 areal_radius = contact_radius |
|
245 end |
|
246 |
|
247 contacts::Array{Int, 1} = zeros(Int, simulation.Nc_max) |
|
248 position_vector = Vector{Vector{Float64}}(undef, simulation.Nc_max) |
|
249 contact_parallel_displacement = |
|
250 Vector{Vector{Float64}}(undef, simulation.Nc_max) |
|
251 contact_rotation = Vector{Vector{Float64}}(undef, simulation.Nc_max) |
|
252 contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max) |
|
253 contact_area::Vector{Float64} = zeros(Float64, simulation.Nc_max) |
|
254 compressive_failure::Vector{Bool} = zeros(Bool, simulation.Nc_max) |
|
255 for i=1:simulation.Nc_max |
|
256 position_vector[i] = zeros(3) |
|
257 contact_rotation[i] = zeros(3) |
|
258 contact_parallel_displacement[i] = zeros(3) |
|
259 end |
|
260 |
|
261 # Create grain object with placeholder values for surface area, volu… |
|
262 # mass, and moment of inertia. |
|
263 grain = GrainCylindrical(density, |
|
264 |
|
265 thickness, |
|
266 contact_radius, |
|
267 areal_radius, |
|
268 1.0, # circumreference |
|
269 1.0, # horizontal_surface_area |
|
270 1.0, # side_surface_area |
|
271 1.0, # volume |
|
272 1.0, # mass |
|
273 1.0, # moment_of_inertia |
|
274 lin_pos, |
|
275 lin_vel, |
|
276 lin_acc, |
|
277 force, |
|
278 [0., 0., 0.], # external_body_force |
|
279 [0., 0., 0.], # lin_disp |
|
280 |
|
281 ang_pos, |
|
282 ang_vel, |
|
283 ang_acc, |
|
284 torque, |
|
285 |
|
286 fixed, |
|
287 allow_x_acc, |
|
288 allow_y_acc, |
|
289 allow_z_acc, |
|
290 rotating, |
|
291 enabled, |
|
292 |
|
293 contact_stiffness_normal, |
|
294 contact_stiffness_tangential, |
|
295 contact_viscosity_normal, |
|
296 contact_viscosity_tangential, |
|
297 contact_static_friction, |
|
298 contact_dynamic_friction, |
|
299 |
|
300 youngs_modulus, |
|
301 poissons_ratio, |
|
302 tensile_strength, |
|
303 shear_strength, |
|
304 strength_heal_rate, |
|
305 fracture_toughness, |
|
306 |
|
307 ocean_drag_coeff_vert, |
|
308 ocean_drag_coeff_horiz, |
|
309 atmosphere_drag_coeff_vert, |
|
310 atmosphere_drag_coeff_horiz, |
|
311 |
|
312 pressure, |
|
313 n_contacts, |
|
314 ocean_grid_pos, |
|
315 atmosphere_grid_pos, |
|
316 contacts, |
|
317 position_vector, |
|
318 contact_parallel_displacement, |
|
319 contact_rotation, |
|
320 contact_age, |
|
321 contact_area, |
|
322 compressive_failure, |
|
323 |
|
324 granular_stress, |
|
325 ocean_stress, |
|
326 atmosphere_stress, |
|
327 |
|
328 thermal_energy, |
|
329 |
|
330 color |
|
331 ) |
|
332 |
|
333 # Overwrite previous placeholder values |
|
334 grain.circumreference = grainCircumreference(grain) |
|
335 grain.horizontal_surface_area = grainHorizontalSurfaceArea(grain) |
|
336 grain.side_surface_area = grainSideSurfaceArea(grain) |
|
337 grain.volume = grainVolume(grain) |
|
338 grain.mass = grainMass(grain) |
|
339 grain.moment_of_inertia = grainMomentOfInertia(grain) |
|
340 |
|
341 # Add to simulation object |
|
342 addGrain!(simulation, grain, verbose) |
|
343 nothing |
|
344 end |
|
345 |
|
346 export grainCircumreference |
|
347 "Returns the circumreference of the grain" |
|
348 function grainCircumreference(grain::GrainCylindrical) |
|
349 return pi*grain.areal_radius*2. |
|
350 end |
|
351 |
|
352 export grainHorizontalSurfaceArea |
|
353 "Returns the top or bottom (horizontal) surface area of the grain" |
|
354 function grainHorizontalSurfaceArea(grain::GrainCylindrical) |
|
355 return pi*grain.areal_radius^2. |
|
356 end |
|
357 |
|
358 export grainSideSurfaceArea |
|
359 "Returns the surface area of the grain sides" |
|
360 function grainSideSurfaceArea(grain::GrainCylindrical) |
|
361 return grainCircumreference(grain)*grain.thickness |
|
362 end |
|
363 |
|
364 export grainVolume |
|
365 "Returns the volume of the grain" |
|
366 function grainVolume(grain::GrainCylindrical) |
|
367 return grainHorizontalSurfaceArea(grain)*grain.thickness |
|
368 end |
|
369 |
|
370 export grainMass |
|
371 "Returns the mass of the grain" |
|
372 function grainMass(grain::GrainCylindrical) |
|
373 return grainVolume(grain)*grain.density |
|
374 end |
|
375 |
|
376 export grainMomentOfInertia |
|
377 "Returns the moment of inertia of the grain" |
|
378 function grainMomentOfInertia(grain::GrainCylindrical) |
|
379 return 0.5*grainMass(grain)*grain.areal_radius^2. |
|
380 end |
|
381 |
|
382 export convertGrainDataToArrays |
|
383 """ |
|
384 Gathers all grain parameters from the `GrainCylindrical` type in continu… |
|
385 arrays in an `GrainArrays` structure. |
|
386 """ |
|
387 function convertGrainDataToArrays(simulation::Simulation) |
|
388 |
|
389 ifarr = GrainArrays( |
|
390 # Material properties |
|
391 ## density |
|
392 Array{Float64}(undef, length(simulation.grains)), |
|
393 |
|
394 # Geometrical properties |
|
395 ## thickness |
|
396 Array{Float64}(undef, length(simulation.grains)), |
|
397 ## contact_radius |
|
398 Array{Float64}(undef, length(simulation.grains)), |
|
399 ## areal_radius |
|
400 Array{Float64}(undef, length(simulation.grains)), |
|
401 ## circumreference |
|
402 Array{Float64}(undef, length(simulation.grains)), |
|
403 ## horizontal_surface_area |
|
404 Array{Float64}(undef, length(simulation.grains)), |
|
405 ## side_surface_area |
|
406 Array{Float64}(undef, length(simulation.grains)), |
|
407 ## volume |
|
408 Array{Float64}(undef, length(simulation.grains)), |
|
409 ## mass |
|
410 Array{Float64}(undef, length(simulation.grains)), |
|
411 ## moment_of_inertia |
|
412 Array{Float64}(undef, length(simulation.grains)), |
|
413 |
|
414 # Linear kinematic degrees of freedom along hori… |
|
415 ## lin_pos |
|
416 zeros(Float64, 3, length(simulation.grains)), |
|
417 ## lin_vel |
|
418 zeros(Float64, 3, length(simulation.grains)), |
|
419 ## lin_acc |
|
420 zeros(Float64, 3, length(simulation.grains)), |
|
421 ## force |
|
422 zeros(Float64, 3, length(simulation.grains)), |
|
423 ## external_body_force |
|
424 zeros(Float64, 3, length(simulation.grains)), |
|
425 ## lin_disp |
|
426 zeros(Float64, 3, length(simulation.grains)), |
|
427 |
|
428 # Angular kinematic degrees of freedom for vert.… |
|
429 ## ang_pos |
|
430 zeros(Float64, 3, length(simulation.grains)), |
|
431 ## ang_vel |
|
432 zeros(Float64, 3, length(simulation.grains)), |
|
433 ## ang_acc |
|
434 zeros(Float64, 3, length(simulation.grains)), |
|
435 ## torque |
|
436 zeros(Float64, 3, length(simulation.grains)), |
|
437 |
|
438 # Kinematic constraint flags |
|
439 ## fixed |
|
440 Array{Int}(undef, length(simulation.grains)), |
|
441 ## allow_x_acc |
|
442 Array{Int}(undef, length(simulation.grains)), |
|
443 ## allow_y_acc |
|
444 Array{Int}(undef, length(simulation.grains)), |
|
445 ## allow_z_acc |
|
446 Array{Int}(undef, length(simulation.grains)), |
|
447 ## rotating |
|
448 Array{Int}(undef, length(simulation.grains)), |
|
449 ## enabled |
|
450 Array{Int}(undef, length(simulation.grains)), |
|
451 |
|
452 # Rheological parameters |
|
453 ## contact_stiffness_normal |
|
454 Array{Float64}(undef, length(simulation.grains)), |
|
455 ## contact_stiffness_tangential |
|
456 Array{Float64}(undef, length(simulation.grains)), |
|
457 ## contact_viscosity_normal |
|
458 Array{Float64}(undef, length(simulation.grains)), |
|
459 ## contact_viscosity_tangential |
|
460 Array{Float64}(undef, length(simulation.grains)), |
|
461 ## contact_static_friction |
|
462 Array{Float64}(undef, length(simulation.grains)), |
|
463 ## contact_dynamic_friction |
|
464 Array{Float64}(undef, length(simulation.grains)), |
|
465 |
|
466 ## youngs_modulus |
|
467 Array{Float64}(undef, length(simulation.grains)), |
|
468 ## poissons_ratio |
|
469 Array{Float64}(undef, length(simulation.grains)), |
|
470 ## tensile_strength |
|
471 Array{Float64}(undef, length(simulation.grains)), |
|
472 ## shear_strength |
|
473 Array{Float64}(undef, length(simulation.grains)), |
|
474 ## strength_heal_rate |
|
475 Array{Float64}(undef, length(simulation.grains)), |
|
476 ## fracture_toughness |
|
477 Array{Float64}(undef, length(simulation.grains)), |
|
478 |
|
479 # Ocean/atmosphere interaction parameters |
|
480 ## ocean_drag_coeff_vert |
|
481 Array{Float64}(undef, length(simulation.grains)), |
|
482 ## ocean_drag_coeff_horiz |
|
483 Array{Float64}(undef, length(simulation.grains)), |
|
484 ## atmosphere_drag_coeff_vert |
|
485 Array{Float64}(undef, length(simulation.grains)), |
|
486 ## atmosphere_drag_coeff_horiz |
|
487 Array{Float64}(undef, length(simulation.grains)), |
|
488 |
|
489 # Interaction |
|
490 ## pressure |
|
491 Array{Float64}(undef, length(simulation.grains)), |
|
492 ## n_contacts |
|
493 Array{Int}(undef, length(simulation.grains)), |
|
494 |
|
495 ## granular_stress |
|
496 zeros(Float64, 3, length(simulation.grains)), |
|
497 ## ocean_stress |
|
498 zeros(Float64, 3, length(simulation.grains)), |
|
499 ## atmosphere_stress |
|
500 zeros(Float64, 3, length(simulation.grains)), |
|
501 |
|
502 ## thermal_energy |
|
503 Array{Float64}(undef, length(simulation.grains)), |
|
504 |
|
505 # Visualization parameters |
|
506 ## color |
|
507 Array{Int}(undef, length(simulation.grains)), |
|
508 |
|
509 ) |
|
510 |
|
511 # fill arrays |
|
512 for i=1:length(simulation.grains) |
|
513 ifarr.density[i] = simulation.grains[i].density |
|
514 |
|
515 ifarr.thickness[i] = simulation.grains[i].thickness |
|
516 ifarr.contact_radius[i] = simulation.grains[i].contact_radius |
|
517 ifarr.areal_radius[i] = simulation.grains[i].areal_radius |
|
518 ifarr.circumreference[i] = simulation.grains[i].circumreference |
|
519 ifarr.horizontal_surface_area[i] = |
|
520 simulation.grains[i].horizontal_surface_area |
|
521 ifarr.side_surface_area[i] = simulation.grains[i].side_surface_a… |
|
522 ifarr.volume[i] = simulation.grains[i].volume |
|
523 ifarr.mass[i] = simulation.grains[i].mass |
|
524 ifarr.moment_of_inertia[i] = simulation.grains[i].moment_of_iner… |
|
525 |
|
526 ifarr.lin_pos[1:3, i] = simulation.grains[i].lin_pos |
|
527 ifarr.lin_vel[1:3, i] = simulation.grains[i].lin_vel |
|
528 ifarr.lin_acc[1:3, i] = simulation.grains[i].lin_acc |
|
529 ifarr.force[1:3, i] = simulation.grains[i].force |
|
530 ifarr.external_body_force[1:3, i] = |
|
531 simulation.grains[i].external_body_force |
|
532 ifarr.lin_disp[1:3, i] = simulation.grains[i].lin_disp |
|
533 |
|
534 ifarr.ang_pos[1:3, i] = simulation.grains[i].ang_pos |
|
535 ifarr.ang_vel[1:3, i] = simulation.grains[i].ang_vel |
|
536 ifarr.ang_acc[1:3, i] = simulation.grains[i].ang_acc |
|
537 ifarr.torque[1:3, i] = simulation.grains[i].torque |
|
538 |
|
539 ifarr.fixed[i] = Int(simulation.grains[i].fixed) |
|
540 ifarr.allow_x_acc[i] = Int(simulation.grains[i].allow_x_acc) |
|
541 ifarr.allow_y_acc[i] = Int(simulation.grains[i].allow_y_acc) |
|
542 ifarr.allow_z_acc[i] = Int(simulation.grains[i].allow_z_acc) |
|
543 ifarr.rotating[i] = Int(simulation.grains[i].rotating) |
|
544 ifarr.enabled[i] = Int(simulation.grains[i].enabled) |
|
545 |
|
546 ifarr.contact_stiffness_normal[i] = |
|
547 simulation.grains[i].contact_stiffness_normal |
|
548 ifarr.contact_stiffness_tangential[i] = |
|
549 simulation.grains[i].contact_stiffness_tangential |
|
550 ifarr.contact_viscosity_normal[i] = |
|
551 simulation.grains[i].contact_viscosity_normal |
|
552 ifarr.contact_viscosity_tangential[i] = |
|
553 simulation.grains[i].contact_viscosity_tangential |
|
554 ifarr.contact_static_friction[i] = |
|
555 simulation.grains[i].contact_static_friction |
|
556 ifarr.contact_dynamic_friction[i] = |
|
557 simulation.grains[i].contact_dynamic_friction |
|
558 |
|
559 ifarr.youngs_modulus[i] = simulation.grains[i].youngs_modulus |
|
560 ifarr.poissons_ratio[i] = simulation.grains[i].poissons_ratio |
|
561 ifarr.tensile_strength[i] = simulation.grains[i].tensile_strength |
|
562 ifarr.shear_strength[i] = simulation.grains[i].shear_strength |
|
563 ifarr.strength_heal_rate[i] = simulation.grains[i].strength_heal… |
|
564 ifarr.fracture_toughness[i] = |
|
565 simulation.grains[i].fracture_toughness |
|
566 |
|
567 ifarr.ocean_drag_coeff_vert[i] = |
|
568 simulation.grains[i].ocean_drag_coeff_vert |
|
569 ifarr.ocean_drag_coeff_horiz[i] = |
|
570 simulation.grains[i].ocean_drag_coeff_horiz |
|
571 ifarr.atmosphere_drag_coeff_vert[i] = |
|
572 simulation.grains[i].atmosphere_drag_coeff_vert |
|
573 ifarr.atmosphere_drag_coeff_horiz[i] = |
|
574 simulation.grains[i].atmosphere_drag_coeff_horiz |
|
575 |
|
576 ifarr.pressure[i] = simulation.grains[i].pressure |
|
577 ifarr.n_contacts[i] = simulation.grains[i].n_contacts |
|
578 |
|
579 ifarr.granular_stress[1:3, i] = simulation.grains[i].granular_st… |
|
580 ifarr.ocean_stress[1:3, i] = simulation.grains[i].ocean_stress |
|
581 ifarr.atmosphere_stress[1:3, i] = simulation.grains[i].atmospher… |
|
582 |
|
583 ifarr.thermal_energy[i] = simulation.grains[i].thermal_energy |
|
584 |
|
585 ifarr.color[i] = simulation.grains[i].color |
|
586 end |
|
587 |
|
588 return ifarr |
|
589 end |
|
590 |
|
591 function deleteGrainArrays!(ifarr::GrainArrays) |
|
592 f1 = zeros(1) |
|
593 f2 = zeros(1,1) |
|
594 i1 = zeros(Int, 1) |
|
595 |
|
596 ifarr.density = f1 |
|
597 |
|
598 ifarr.thickness = f1 |
|
599 ifarr.contact_radius = f1 |
|
600 ifarr.areal_radius = f1 |
|
601 ifarr.circumreference = f1 |
|
602 ifarr.horizontal_surface_area = f1 |
|
603 ifarr.side_surface_area = f1 |
|
604 ifarr.volume = f1 |
|
605 ifarr.mass = f1 |
|
606 ifarr.moment_of_inertia = f1 |
|
607 |
|
608 ifarr.lin_pos = f2 |
|
609 ifarr.lin_vel = f2 |
|
610 ifarr.lin_acc = f2 |
|
611 ifarr.force = f2 |
|
612 ifarr.external_body_force = f2 |
|
613 ifarr.lin_disp = f2 |
|
614 |
|
615 ifarr.ang_pos = f2 |
|
616 ifarr.ang_vel = f2 |
|
617 ifarr.ang_acc = f2 |
|
618 ifarr.torque = f2 |
|
619 |
|
620 ifarr.fixed = i1 |
|
621 ifarr.allow_x_acc = i1 |
|
622 ifarr.allow_y_acc = i1 |
|
623 ifarr.allow_z_acc = i1 |
|
624 ifarr.rotating = i1 |
|
625 ifarr.enabled = i1 |
|
626 |
|
627 ifarr.contact_stiffness_normal = f1 |
|
628 ifarr.contact_stiffness_tangential = f1 |
|
629 ifarr.contact_viscosity_normal = f1 |
|
630 ifarr.contact_viscosity_tangential = f1 |
|
631 ifarr.contact_static_friction = f1 |
|
632 ifarr.contact_dynamic_friction = f1 |
|
633 |
|
634 ifarr.youngs_modulus = f1 |
|
635 ifarr.poissons_ratio = f1 |
|
636 ifarr.tensile_strength = f1 |
|
637 ifarr.shear_strength = f1 |
|
638 ifarr.strength_heal_rate = f1 |
|
639 ifarr.fracture_toughness = f1 |
|
640 |
|
641 ifarr.ocean_drag_coeff_vert = f1 |
|
642 ifarr.ocean_drag_coeff_horiz = f1 |
|
643 ifarr.atmosphere_drag_coeff_vert = f1 |
|
644 ifarr.atmosphere_drag_coeff_horiz = f1 |
|
645 |
|
646 ifarr.pressure = f1 |
|
647 ifarr.n_contacts = i1 |
|
648 |
|
649 ifarr.granular_stress = f2 |
|
650 ifarr.ocean_stress = f2 |
|
651 ifarr.atmosphere_stress = f2 |
|
652 |
|
653 ifarr.thermal_energy = f1 |
|
654 |
|
655 ifarr.color = i1 |
|
656 nothing |
|
657 end |
|
658 |
|
659 export printGrainInfo |
|
660 """ |
|
661 printGrainInfo(grain::GrainCylindrical) |
|
662 |
|
663 Prints the contents of an grain to stdout in a formatted style. |
|
664 """ |
|
665 function printGrainInfo(f::GrainCylindrical) |
|
666 println(" density: $(f.density) kg/m^3") |
|
667 println(" thickness: $(f.thickness) m") |
|
668 println(" contact_radius: $(f.contact_radius) m") |
|
669 println(" areal_radius: $(f.areal_radius) m") |
|
670 println(" circumreference: $(f.circumreference) m") |
|
671 println(" horizontal_surface_area: $(f.horizontal_surface_area) m^2… |
|
672 println(" side_surface_area: $(f.side_surface_area) m^2") |
|
673 println(" volume: $(f.volume) m^3") |
|
674 println(" mass: $(f.mass) kg") |
|
675 println(" moment_of_inertia: $(f.moment_of_inertia) kg*m^2\n") |
|
676 |
|
677 println(" lin_pos: $(f.lin_pos) m") |
|
678 println(" lin_vel: $(f.lin_vel) m/s") |
|
679 println(" lin_acc: $(f.lin_acc) m/s^2") |
|
680 println(" force: $(f.force) N\n") |
|
681 println(" external_body_force: $(f.external_body_force) N\n") |
|
682 |
|
683 println(" ang_pos: $(f.ang_pos) rad") |
|
684 println(" ang_vel: $(f.ang_vel) rad/s") |
|
685 println(" ang_acc: $(f.ang_acc) rad/s^2") |
|
686 println(" torque: $(f.torque) N*m\n") |
|
687 |
|
688 println(" fixed: $(f.fixed)") |
|
689 println(" allow_x_acc: $(f.allow_x_acc)") |
|
690 println(" allow_y_acc: $(f.allow_y_acc)") |
|
691 println(" allow_z_acc: $(f.allow_z_acc)") |
|
692 println(" rotating: $(f.rotating)") |
|
693 println(" enabled: $(f.enabled)\n") |
|
694 |
|
695 println(" k_n: $(f.contact_stiffness_normal) N/m") |
|
696 println(" k_t: $(f.contact_stiffness_tangential) N/m") |
|
697 println(" γ_n: $(f.contact_viscosity_normal) N/(m/s)") |
|
698 println(" γ_t: $(f.contact_viscosity_tangential) N/(m/s)") |
|
699 println(" μ_s: $(f.contact_static_friction)") |
|
700 println(" μ_d: $(f.contact_dynamic_friction)\n") |
|
701 |
|
702 println(" E: $(f.youngs_modulus) Pa") |
|
703 println(" ν: $(f.poissons_ratio)") |
|
704 println(" tensile_strength: $(f.tensile_strength) Pa") |
|
705 println(" shear_strength: $(f.shear_strength) Pa") |
|
706 println(" strength_heal_rate: $(f.strength_heal_rate) Pa/s") |
|
707 println(" fracture_toughness: $(f.fracture_toughness) m^0.5 Pa\n") |
|
708 |
|
709 println(" c_o_v: $(f.ocean_drag_coeff_vert)") |
|
710 println(" c_o_h: $(f.ocean_drag_coeff_horiz)") |
|
711 println(" c_a_v: $(f.atmosphere_drag_coeff_vert)") |
|
712 println(" c_a_h: $(f.atmosphere_drag_coeff_horiz)\n") |
|
713 |
|
714 println(" pressure: $(f.pressure) Pa") |
|
715 println(" n_contacts: $(f.n_contacts)") |
|
716 println(" contacts: $(f.contacts)") |
|
717 println(" pos_ij: $(f.position_vector)\n") |
|
718 println(" δ_t: $(f.contact_parallel_displacement)\n") |
|
719 println(" θ_t: $(f.contact_rotation)\n") |
|
720 |
|
721 println(" granular_stress: $(f.granular_stress) Pa") |
|
722 println(" ocean_stress: $(f.ocean_stress) Pa") |
|
723 println(" atmosphere_stress: $(f.atmosphere_stress) Pa\n") |
|
724 |
|
725 println(" thermal_energy: $(f.thermal_energy) J\n") |
|
726 |
|
727 println(" color: $(f.color)\n") |
|
728 nothing |
|
729 end |
|
730 |
|
731 export grainKineticTranslationalEnergy |
|
732 "Returns the translational kinetic energy of the grain" |
|
733 function grainKineticTranslationalEnergy(grain::GrainCylindrical) |
|
734 return 0.5*grain.mass*norm(grain.lin_vel)^2. |
|
735 end |
|
736 |
|
737 export totalGrainKineticTranslationalEnergy |
|
738 """ |
|
739 totalGrainKineticTranslationalEnergy(simulation) |
|
740 |
|
741 Returns the sum of translational kinetic energies of all grains in a |
|
742 simulation |
|
743 """ |
|
744 function totalGrainKineticTranslationalEnergy(simulation::Simulation) |
|
745 E_sum = 0. |
|
746 for grain in simulation.grains |
|
747 E_sum += grainKineticTranslationalEnergy(grain) |
|
748 end |
|
749 return E_sum |
|
750 end |
|
751 |
|
752 export grainKineticRotationalEnergy |
|
753 "Returns the rotational kinetic energy of the grain" |
|
754 function grainKineticRotationalEnergy(grain::GrainCylindrical) |
|
755 return 0.5*grain.moment_of_inertia*norm(grain.ang_vel)^2. |
|
756 end |
|
757 |
|
758 export totalGrainKineticRotationalEnergy |
|
759 """ |
|
760 totalGrainKineticRotationalEnergy(simulation) |
|
761 |
|
762 Returns the sum of rotational kinetic energies of all grains in a simula… |
|
763 """ |
|
764 function totalGrainKineticRotationalEnergy(simulation::Simulation) |
|
765 E_sum = 0. |
|
766 for grain in simulation.grains |
|
767 E_sum += grainKineticRotationalEnergy(grain) |
|
768 end |
|
769 return E_sum |
|
770 end |
|
771 |
|
772 export grainThermalEnergy |
|
773 "Returns the thermal energy of the grain, produced by Coulomb slip" |
|
774 function grainThermalEnergy(grain::GrainCylindrical) |
|
775 return grain.thermal_energy |
|
776 end |
|
777 |
|
778 export totalGrainThermalEnergy |
|
779 """ |
|
780 totalGrainKineticTranslationalEnergy(simulation) |
|
781 |
|
782 Returns the sum of thermal energy of all grains in a simulation |
|
783 """ |
|
784 function totalGrainThermalEnergy(simulation::Simulation) |
|
785 E_sum = 0. |
|
786 for grain in simulation.grains |
|
787 E_sum += grainThermalEnergy(grain) |
|
788 end |
|
789 return E_sum |
|
790 end |
|
791 |
|
792 export addBodyForce! |
|
793 """ |
|
794 setBodyForce!(grain, force) |
|
795 |
|
796 Add to the value of the external body force on a grain. |
|
797 |
|
798 # Arguments |
|
799 * `grain::GrainCylindrical`: the grain to set the body force for. |
|
800 * `force::Vector{Float64}`: a vector of force [N] |
|
801 """ |
|
802 function addBodyForce!(grain::GrainCylindrical, force::Vector{Float64}) |
|
803 grain.external_body_force += vecTo3d(force) |
|
804 end |
|
805 |
|
806 export setBodyForce! |
|
807 """ |
|
808 setBodyForce!(grain, force) |
|
809 |
|
810 Set the value of the external body force on a grain. |
|
811 |
|
812 # Arguments |
|
813 * `grain::GrainCylindrical`: the grain to set the body force for. |
|
814 * `force::Vector{Float64}`: a vector of force [N] |
|
815 """ |
|
816 function setBodyForce!(grain::GrainCylindrical, force::Vector{Float64}) |
|
817 grain.external_body_force = vecTo3d(force) |
|
818 end |
|
819 |
|
820 export compareGrains |
|
821 """ |
|
822 compareGrains(if1::GrainCylindrical, if2::GrainCylindrical) |
|
823 |
|
824 Compare values of two grain objects using the `Base.Test` framework. |
|
825 """ |
|
826 function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical) |
|
827 |
|
828 @test if1.density ≈ if2.density |
|
829 @test if1.thickness ≈ if2.thickness |
|
830 @test if1.contact_radius ≈ if2.contact_radius |
|
831 @test if1.areal_radius ≈ if2.areal_radius |
|
832 @test if1.circumreference ≈ if2.circumreference |
|
833 @test if1.horizontal_surface_area ≈ if2.horizontal_surface_area |
|
834 @test if1.side_surface_area ≈ if2.side_surface_area |
|
835 @test if1.volume ≈ if2.volume |
|
836 @test if1.mass ≈ if2.mass |
|
837 @test if1.moment_of_inertia ≈ if2.moment_of_inertia |
|
838 |
|
839 @test if1.lin_pos ≈ if2.lin_pos |
|
840 @test if1.lin_vel ≈ if2.lin_vel |
|
841 @test if1.lin_acc ≈ if2.lin_acc |
|
842 @test if1.force ≈ if2.force |
|
843 @test if1.external_body_force ≈ if2.external_body_force |
|
844 @test if1.lin_disp ≈ if2.lin_disp |
|
845 |
|
846 @test if1.ang_pos ≈ if2.ang_pos |
|
847 @test if1.ang_vel ≈ if2.ang_vel |
|
848 @test if1.ang_acc ≈ if2.ang_acc |
|
849 @test if1.torque ≈ if2.torque |
|
850 |
|
851 @test if1.fixed == if2.fixed |
|
852 @test if1.rotating == if2.rotating |
|
853 @test if1.enabled == if2.enabled |
|
854 |
|
855 @test if1.contact_stiffness_normal ≈ if2.contact_stiffness_normal |
|
856 @test if1.contact_stiffness_tangential ≈ if2.contact_stiffness_tan… |
|
857 @test if1.contact_viscosity_normal ≈ if2.contact_viscosity_normal |
|
858 @test if1.contact_viscosity_tangential ≈ if2.contact_viscosity_tan… |
|
859 @test if1.contact_static_friction ≈ if2.contact_static_friction |
|
860 @test if1.contact_dynamic_friction ≈ if2.contact_dynamic_friction |
|
861 |
|
862 @test if1.youngs_modulus ≈ if2.youngs_modulus |
|
863 @test if1.poissons_ratio ≈ if2.poissons_ratio |
|
864 @test if1.tensile_strength ≈ if2.tensile_strength |
|
865 @test if1.shear_strength ≈ if2.shear_strength |
|
866 @test if1.strength_heal_rate ≈ if2.strength_heal_rate |
|
867 @test if1.fracture_toughness ≈ if2.fracture_toughness |
|
868 |
|
869 @test if1.ocean_drag_coeff_vert ≈ if2.ocean_drag_coeff_vert |
|
870 @test if1.ocean_drag_coeff_horiz ≈ if2.ocean_drag_coeff_horiz |
|
871 @test if1.atmosphere_drag_coeff_vert ≈ if2.atmosphere_drag_coeff_v… |
|
872 @test if1.atmosphere_drag_coeff_horiz ≈ if2.atmosphere_drag_coeff_… |
|
873 |
|
874 @test if1.pressure ≈ if2.pressure |
|
875 @test if1.n_contacts == if2.n_contacts |
|
876 @test if1.ocean_grid_pos == if2.ocean_grid_pos |
|
877 @test if1.contacts == if2.contacts |
|
878 @test if1.position_vector == if2.position_vector |
|
879 @test if1.contact_parallel_displacement == if2.contact_parallel_disp… |
|
880 @test if1.contact_rotation == if2.contact_rotation |
|
881 @test if1.contact_age ≈ if2.contact_age |
|
882 @test if1.contact_area ≈ if2.contact_area |
|
883 @test if1.compressive_failure ≈ if2.compressive_failure |
|
884 |
|
885 @test if1.granular_stress ≈ if2.granular_stress |
|
886 @test if1.ocean_stress ≈ if2.ocean_stress |
|
887 @test if1.atmosphere_stress ≈ if2.atmosphere_stress |
|
888 |
|
889 @test if1.thermal_energy ≈ if2.thermal_energy |
|
890 |
|
891 @test if1.color ≈ if2.color |
|
892 nothing |
|
893 end |
|
894 |
|
895 export enableOceanDrag! |
|
896 """ |
|
897 enableOceanDrag!(grain) |
|
898 |
|
899 Enable ocean-caused drag on the grain, with values by Hunke and Comeau (… |
|
900 """ |
|
901 function enableOceanDrag!(grain::GrainCylindrical) |
|
902 grain.ocean_drag_coeff_vert = 0.85 |
|
903 grain.ocean_drag_coeff_horiz = 5e-4 |
|
904 end |
|
905 |
|
906 export enableAtmosphereDrag! |
|
907 """ |
|
908 enableAtmosphereDrag!(grain) |
|
909 |
|
910 Enable atmosphere-caused drag on the grain, with values by Hunke and Com… |
|
911 (2011). |
|
912 """ |
|
913 function enableAtmosphereDrag!(grain::GrainCylindrical) |
|
914 grain.atmosphere_drag_coeff_vert = 0.4 |
|
915 grain.atmosphere_drag_coeff_horiz = 2.5e-4 |
|
916 end |
|
917 export disableOceanDrag! |
|
918 """ |
|
919 disableOceanDrag!(grain) |
|
920 |
|
921 Disable ocean-caused drag on the grain. |
|
922 """ |
|
923 function disableOceanDrag!(grain::GrainCylindrical) |
|
924 grain.ocean_drag_coeff_vert = 0. |
|
925 grain.ocean_drag_coeff_horiz = 0. |
|
926 end |
|
927 |
|
928 export disableAtmosphereDrag! |
|
929 """ |
|
930 disableAtmosphereDrag!(grain) |
|
931 |
|
932 Disable atmosphere-caused drag on the grain. |
|
933 """ |
|
934 function disableAtmosphereDrag!(grain::GrainCylindrical) |
|
935 grain.atmosphere_drag_coeff_vert = 0. |
|
936 grain.atmosphere_drag_coeff_horiz = 0. |
|
937 end |
|
938 |
|
939 export zeroKinematics! |
|
940 """ |
|
941 zeroKinematics!(simulation) |
|
942 |
|
943 Set all grain forces, torques, accelerations, and velocities (linear and |
|
944 rotational) to zero in order to get rid of all kinetic energy. |
|
945 """ |
|
946 function zeroKinematics!(sim::Simulation) |
|
947 for grain in sim.grains |
|
948 grain.lin_vel .= zeros(3) |
|
949 grain.lin_acc .= zeros(3) |
|
950 grain.force .= zeros(3) |
|
951 grain.lin_disp .= zeros(3) |
|
952 grain.ang_vel .= zeros(3) |
|
953 grain.ang_acc .= zeros(3) |
|
954 grain.torque .= zeros(3) |
|
955 end |
|
956 end |
