grain.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