tgrain.jl - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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, lin_acc,
10 force, ang_pos, ang_vel, ang_acc, torque,
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_acc,
27 rotating, enabled, verbose,
28 ocean_grid_pos, atmosphere_grid_pos,
29 n_contact, granular_stress, ocean_stress,
30 atmosphere_stress,
31 thermal_energy,
32 color])
33
34 Creates and adds a cylindrical grain to a simulation. Most of the arguments
35 are optional, and come with default values. The only required arguments are
36 `simulation`, `lin_pos`, `contact_radius`, and `thickness`.
37
38 # Arguments
39 * `simulation::Simulation`: the simulation object where the grain should be
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 *y*, 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 concentration
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 *y*, and
50 the *z* component will be set to zero.
51 * `lin_acc::Vector{Float64} = [0., 0., 0.]`: linear acceleration [m/s^2]. If a
52 two-component vector is used, the values will be mapped to *x* and *y*, 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 *y*, and
56 the *z* component will be set to zero.
57 * `ang_pos::Float64 = [0., 0., 0.]`: angular position around its center vertical
58 axis [rad]. If a scalar is used, the value will be mapped to *z*, and the
59 *x* and *y* components will be set to zero.
60 * `ang_vel::Float64 = [0., 0., 0.]`: angular velocity around its center vertical
61 axis [rad/s]. If a scalar is used, the value will be mapped to *z*, and the
62 *x* and *y* components will be set to zero.
63 * `ang_acc::Float64 = [0., 0., 0.]`: angular acceleration around its center
64 vertical axis [rad/s^2]. If a scalar is used, the value will be mapped to
65 *z*, and the *x* and *y* components will be set to zero.
66 * `torque::Float64 = [0., 0., 0.]`: torque around its center vertical axis
67 [N*m]. If a scalar is used, the value will be mapped to *z*, and the *x* and
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/m];
71 overridden if `youngs_modulus` is set to a positive value.
72 * `contact_stiffness_tangential::Float64 = 0.`: contact-tangential stiffness
73 [N/m]; overridden if `youngs_modulus` is set to a positive value.
74 * `contact_viscosity_normal::Float64 = 0.`: contact-normal viscosity [N/m/s].
75 * `contact_viscosity_tangential::Float64 = 0.`: contact-tangential viscosity
76 [N/m/s].
77 * `contact_static_friction::Float64 = 0.4`: contact static Coulomb frictional
78 coefficient [-].
79 * `contact_dynamic_friction::Float64 = 0.4`: contact dynamic Coulomb frictional
80 coefficient [-].
81 * `youngs_modulus::Float64 = 2e7`: elastic modulus [Pa]; overrides any value
82 set for `contact_stiffness_normal`.
83 * `poissons_ratio::Float64 = 0.185`: Poisson's ratio, used to determine the
84 contact-tangential stiffness from `youngs_modulus` [-].
85 * `tensile_strength::Float64 = 0.`: contact-tensile (cohesive) bond strength
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 influences the
91 maximum compressive strength on granular contact [m^{1/2}*Pa]. A value
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 ocean
94 against grain sides [-].
95 * `ocean_drag_coeff_horiz::Float64 = 5e-4`: horizontal drag coefficient for
96 ocean against grain bottom [-].
97 * `atmosphere_drag_coeff_vert::Float64 = 0.4`: vertical drag coefficient for
98 atmosphere against grain sides [-].
99 * `atmosphere_drag_coeff_horiz::Float64 = 2.5e-4`: horizontal drag coefficient
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. zero).
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 function
109 call.
110 * `ocean_grid_pos::Array{Int, 1} = [0, 0]`: position of grain in the ocean
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 grain
116 from granular interactions [Pa].
117 * `ocean_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on grain from
118 ocean drag [Pa].
119 * `atmosphere_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on grain
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 the grain
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.5`
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 passed.
134
135 The following example will create a grain with tensile and shear strength, and a
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 useful
146 for creating regular arrangements:
147
148 ```julia
149 for i=1:5
150 Granular.addGrainCylindrical!(sim, [i*2., 0., 3.], 1., .5, fixed=true)
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 = 0.,
169 contact_viscosity_normal::Float64 = 0.,
170 contact_viscosity_tangential::Float64 = 0.,
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, # Hopkins 2004
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, # H&C 2011
180 ocean_drag_coeff_horiz::Float64 = 5e-4, # H&C 2011
181 atmosphere_drag_coeff_vert::Float64 = 0.4, # H&C 2011
182 atmosphere_drag_coeff_horiz::Float64 = 2.5e-4, # H&C2011
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, 0],
193 n_contacts::Int = 0,
194 granular_stress::Vector{Float64} = [0., 0., 0.],
195 ocean_stress::Vector{Float64} = [0., 0., 0.],
196 atmosphere_stress::Vector{Float64} = [0., 0., 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, volume,
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 continuous
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 horiz plane
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. rot.
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_area
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_inertia
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_rate
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_stress
580 ifarr.ocean_stress[1:3, i] = simulation.grains[i].ocean_stress
581 ifarr.atmosphere_stress[1:3, i] = simulation.grains[i].atmosphere_stress
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 simulation
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_tangential
857 @test if1.contact_viscosity_normal ≈ if2.contact_viscosity_normal
858 @test if1.contact_viscosity_tangential ≈ if2.contact_viscosity_tangential
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_vert
872 @test if1.atmosphere_drag_coeff_horiz ≈ if2.atmosphere_drag_coeff_horiz
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_displacement
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 (2011).
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 Comeau
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