interaction.jl (15104B)
1 ## Interaction functions 2 using LinearAlgebra 3 4 export interact! 5 """ 6 interact!(simulation::Simulation) 7 8 Resolve mechanical interaction between all particle pairs. 9 10 # Arguments 11 * `simulation::Simulation`: the simulation object containing the grains. 12 """ 13 function interact!(simulation::Simulation) 14 for i=1:length(simulation.grains) 15 for ic=1:simulation.Nc_max 16 17 j = simulation.grains[i].contacts[ic] 18 19 if i > j # skip i > j and j == 0 20 continue 21 end 22 23 interactGrains!(simulation, i, j, ic) 24 end 25 end 26 27 for i=1:length(simulation.grains) 28 @inbounds simulation.grains[i].granular_stress = 29 simulation.grains[i].force/ 30 simulation.grains[i].horizontal_surface_area 31 end 32 nothing 33 end 34 35 36 """ 37 interactWalls!(sim) 38 39 Find and resolve interactions between the dynamic walls (`simulation.walls`) and 40 the grains. The contact model uses linear elasticity, with stiffness based on 41 the grain Young's modulus `grian.E` or elastic stiffness `grain.k_n`. The 42 interaction is frictionless in the tangential direction. 43 44 # Arguments 45 * `simulation::Simulation`: the simulation object containing the grains and 46 dynamic walls. 47 """ 48 function interactWalls!(sim::Simulation) 49 50 orientation::Float64 = 0.0 51 δ_n::Float64 = 0.0 52 k_n::Float64 = 0.0 53 γ_n::Float64 = 0.0 54 vel_n::Float64 = 0.0 55 force_n::Float64 = 0.0 56 57 for iw=1:length(sim.walls) 58 for i=1:length(sim.grains) 59 60 orientation = sign(dot(sim.walls[iw].normal, 61 sim.grains[i].lin_pos) - sim.walls[iw].pos) 62 63 # get overlap distance by projecting grain position onto wall-normal 64 # vector. Overlap when δ_n < 0.0 65 δ_n = norm(dot(sim.walls[iw].normal[1:2], 66 sim.grains[i].lin_pos[1:2]) - 67 sim.walls[iw].pos) - sim.grains[i].contact_radius 68 69 vel_n = dot(sim.walls[iw].normal[1:2], sim.grains[i].lin_vel[1:2]) 70 71 if δ_n < 0. 72 if sim.grains[i].youngs_modulus > 0. 73 k_n = sim.grains[i].youngs_modulus * sim.grains[i].thickness 74 else 75 k_n = sim.grains[i].contact_stiffness_normal 76 end 77 78 γ_n = sim.walls[iw].contact_viscosity_normal 79 80 if k_n ≈ 0. && γ_n ≈ 0. # No interaction 81 force_n = 0. 82 83 elseif k_n > 0. && γ_n ≈ 0. # Elastic (Hookean) 84 force_n = k_n * δ_n 85 86 elseif k_n > 0. && γ_n > 0. # Elastic-viscous (Kelvin-Voigt) 87 force_n = k_n * δ_n + γ_n * vel_n 88 89 # Make sure that the visous component doesn't dominate over 90 # elasticity 91 if force_n > 0. 92 force_n = 0. 93 end 94 95 else 96 error("unknown contact_normal_rheology " * 97 "(k_n = $k_n, γ_n = $γ_n") 98 end 99 100 sim.grains[i].force += -force_n .* sim.walls[iw].normal .* 101 orientation 102 sim.walls[iw].force += force_n * orientation 103 end 104 end 105 end 106 end 107 108 export interactGrains! 109 """ 110 interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int) 111 112 Resolve an grain-to-grain interaction using a prescibed contact law. This 113 function adds the compressive force of the interaction to the grain 114 `pressure` field of mean compressive stress on the grain sides. 115 116 The function uses the macroscopic contact-stiffness parameterization based on 117 Young's modulus and Poisson's ratio if Young's modulus is a positive value. 118 """ 119 function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int) 120 121 if !simulation.grains[i].enabled || !simulation.grains[j].enabled 122 return 123 end 124 125 force_n = 0. # Contact-normal force 126 force_t = 0. # Contact-parallel (tangential) force 127 128 # Inter-position vector, use stored value which is corrected for boundary 129 # periodicity 130 p = simulation.grains[i].position_vector[ic][1:2] 131 dist = norm(p) 132 133 r_i = simulation.grains[i].contact_radius 134 r_j = simulation.grains[j].contact_radius 135 136 # Floe distance; <0: compression, >0: tension 137 δ_n = dist - (r_i + r_j) 138 139 # Local axes 140 n = p/max(dist, 1e-12) 141 t = [-n[2], n[1]] 142 143 # Contact kinematics (2d) 144 vel_lin = simulation.grains[i].lin_vel[1:2] - 145 simulation.grains[j].lin_vel[1:2] 146 vel_n = dot(vel_lin, n) 147 148 vel_t = dot(vel_lin, t) - 149 harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel[3] + 150 simulation.grains[j].ang_vel[3]) 151 152 # Correct old tangential displacement for contact rotation and add new 153 δ_t0 = simulation.grains[i].contact_parallel_displacement[ic][1:2] 154 if !simulation.grains[i].compressive_failure[ic] 155 δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step 156 else 157 # Use δ_t as total slip-distance vector 158 δ_t = δ_t0 .+ (vel_lin .+ vel_t.*t) .* simulation.time_step 159 end 160 161 # Determine the contact rotation (2d) 162 θ_t = simulation.grains[i].contact_rotation[ic][3] + 163 (simulation.grains[j].ang_vel[3] - 164 simulation.grains[i].ang_vel[3]) * simulation.time_step 165 166 # Effective radius 167 R_ij = harmonicMean(r_i, r_j) 168 169 # Determine which ice floe is the thinnest 170 h_min, idx_thinnest = findmin([simulation.grains[i].thickness, 171 simulation.grains[j].thickness]) 172 idx_thinnest == 1 ? idx_thinnest = i : idx_thinnest = j 173 174 # Contact length along z 175 Lz_ij = h_min 176 177 # Contact area 178 A_ij = R_ij*Lz_ij 179 simulation.grains[i].contact_area[ic] = A_ij 180 181 # Contact mechanical parameters 182 if simulation.grains[i].youngs_modulus > 0. && 183 simulation.grains[j].youngs_modulus > 0. 184 185 E = harmonicMean(simulation.grains[i].youngs_modulus, 186 simulation.grains[j].youngs_modulus) 187 ν = harmonicMean(simulation.grains[i].poissons_ratio, 188 simulation.grains[j].poissons_ratio) 189 190 # Effective normal and tangential stiffness 191 k_n = E * A_ij/R_ij 192 #k_t = k_n*ν # Kneib et al 2016 193 k_t = k_n * 2. * (1. - ν^2.) / ((2. - ν) * (1. + ν)) # Obermayr 2011 194 195 else # Micromechanical parameterization: k_n and k_t set explicitly 196 k_n = harmonicMean(simulation.grains[i].contact_stiffness_normal, 197 simulation.grains[j].contact_stiffness_normal) 198 199 k_t = harmonicMean(simulation.grains[i].contact_stiffness_tangential, 200 simulation.grains[j].contact_stiffness_tangential) 201 end 202 203 γ_n = harmonicMean(simulation.grains[i].contact_viscosity_normal, 204 simulation.grains[j].contact_viscosity_normal) 205 206 γ_t = harmonicMean(simulation.grains[i].contact_viscosity_tangential, 207 simulation.grains[j].contact_viscosity_tangential) 208 209 μ_d_minimum = min(simulation.grains[i].contact_dynamic_friction, 210 simulation.grains[j].contact_dynamic_friction) 211 212 # Determine contact forces 213 if k_n ≈ 0. && γ_n ≈ 0. # No interaction 214 force_n = 0. 215 216 elseif k_n > 0. && γ_n ≈ 0. # Elastic (Hookean) 217 force_n = -k_n*δ_n 218 219 elseif k_n > 0. && γ_n > 0. # Elastic-viscous (Kelvin-Voigt) 220 force_n = -k_n*δ_n - γ_n*vel_n 221 if force_n < 0. 222 force_n = 0. 223 end 224 225 else 226 error("unknown contact_normal_rheology (k_n = $k_n, γ_n = $γ_n, 227 E = $E, A_ij = $A_ij, R_ij = $R_ij) 228 ") 229 end 230 231 # Determine the compressive strength in Pa by the contact thickness and the 232 # minimum compressive strength [N] 233 compressive_strength = min(simulation.grains[i].fracture_toughness, 234 simulation.grains[j].fracture_toughness) * 235 Lz_ij^1.5 236 237 # Grain-pair moment of inertia [m^4] 238 I_ij = 2.0/3.0*R_ij^3*min(simulation.grains[i].thickness, 239 simulation.grains[j].thickness) 240 241 # Contact tensile strength increases linearly with contact age until 242 # tensile stress exceeds tensile strength. 243 tensile_strength = min(simulation.grains[i].contact_age[ic]* 244 simulation.grains[i].strength_heal_rate, 245 simulation.grains[i].tensile_strength) 246 shear_strength = min(simulation.grains[i].contact_age[ic]* 247 simulation.grains[i].strength_heal_rate, 248 simulation.grains[i].shear_strength) 249 M_t = 0.0 250 if tensile_strength > 0.0 251 # Determine bending momentum on contact [N*m], 252 # (converting k_n to E to bar(k_n)) 253 M_t = (k_n*R_ij/(A_ij*(simulation.grains[i].contact_radius + 254 simulation.grains[j].contact_radius)))*I_ij*θ_t 255 end 256 257 # Reset contact age (breaking bond) if bond strength is exceeded 258 259 if δ_n >= 0.0 && norm(force_n)/A_ij + norm(M_t)*R_ij/I_ij > tensile_strength 260 force_n = 0. 261 force_t = 0. 262 simulation.grains[i].contacts[ic] = 0 # remove contact 263 simulation.grains[i].n_contacts -= 1 264 simulation.grains[j].n_contacts -= 1 265 end 266 267 if !simulation.grains[i].compressive_failure[ic] 268 if k_t ≈ 0. && γ_t ≈ 0. 269 # do nothing 270 271 elseif k_t ≈ 0. && γ_t > 0. 272 force_t = norm(γ_t * vel_t) 273 274 # Coulomb slip 275 if force_t > μ_d_minimum*norm(force_n) 276 force_t = μ_d_minimum*norm(force_n) 277 278 # Nguyen et al 2009 "Thermomechanical modelling of friction effects 279 # in granular flows using the discrete element method" 280 E_shear = norm(force_t)*norm(vel_t)*simulation.time_step 281 282 # Assume equal thermal properties 283 simulation.grains[i].thermal_energy += 0.5*E_shear 284 simulation.grains[j].thermal_energy += 0.5*E_shear 285 end 286 if vel_t > 0. 287 force_t = -force_t 288 end 289 290 elseif k_t > 0. 291 292 force_t = -k_t*δ_t - γ_t*vel_t 293 294 # Coulomb slip 295 if norm(force_t) > μ_d_minimum*norm(force_n) 296 force_t = μ_d_minimum*norm(force_n)*force_t/norm(force_t) 297 δ_t = (-force_t - γ_t*vel_t)/k_t 298 299 # Nguyen et al 2009 "Thermomechanical modelling of friction 300 # effects in granular flows using the discrete element method" 301 E_shear = norm(force_t)*norm(vel_t)*simulation.time_step 302 303 # Assume equal thermal properties 304 simulation.grains[i].thermal_energy += 0.5*E_shear 305 simulation.grains[j].thermal_energy += 0.5*E_shear 306 end 307 308 else 309 error("unknown contact_tangential_rheology (k_t = $k_t, γ_t = $γ_t") 310 end 311 end 312 313 just_failed = false 314 # Detect compressive failure if compressive force exceeds compressive 315 # strength 316 if δ_n <= 0.0 && compressive_strength > 0. && 317 !simulation.grains[i].compressive_failure[ic] && 318 norm(force_n.*n .+ force_t.*t) >= compressive_strength 319 320 # Register that compressive failure has occurred for this contact 321 simulation.grains[i].compressive_failure[ic] = true 322 just_failed = true 323 end 324 325 # Overwrite previous force results if compressive failure occurred 326 if simulation.grains[i].compressive_failure[ic] && δ_n < 0.0 327 328 # Normal stress on horizontal interface, assuming bending stresses are 329 # negligible (ocean density and gravity are hardcoded for now) 330 σ_n = (1e3 - simulation.grains[i].density) * 9.81 * 331 (simulation.grains[i].thickness + simulation.grains[j].thickness) 332 333 # Horizontal overlap area (two overlapping circles) 334 if 2.0*min(r_i, r_j) <= abs(δ_n) 335 336 # Complete overlap 337 A_h = π*min(r_i, r_j)^2 338 339 else 340 # Partial overlap, eq 14 from 341 # http://mathworld.wolfram.com/Circle-CircleIntersection.html 342 A_h = r_i^2.0*acos((dist^2.0 + r_i^2.0 - r_j^2.0)/(2.0*dist*r_i)) + 343 r_j^2.0*acos((dist^2.0 + r_j^2.0 - r_i^2.0)/(2.0*dist*r_j)) - 344 0.5*sqrt((-dist + r_i + r_j)* 345 ( dist + r_i - r_j)* 346 ( dist - r_i + r_j)* 347 ( dist + r_i + r_j)) 348 end 349 simulation.grains[i].contact_area[ic] = A_h 350 351 if just_failed 352 # Use δ_t as travel distance on horizontal slip surface, and set the 353 # original displacement to Coulomb-frictional limit in the direction 354 # of the pre-failure force 355 F = force_n.*n .+ force_t.*t 356 σ_t = μ_d_minimum*norm(σ_n).*F/norm(F) 357 δ_t = σ_t*A_h/(-k_t) 358 359 # Save energy loss in E_shear 360 E_shear = norm(F)*norm(vel_lin)*simulation.time_step 361 simulation.grains[i].thermal_energy += 0.5*E_shear 362 simulation.grains[j].thermal_energy += 0.5*E_shear 363 364 else 365 366 # Find horizontal stress on slip interface 367 σ_t = -k_t*δ_t/A_h 368 369 # Check for Coulomb-failure on the slip interface 370 if norm(σ_t) > μ_d_minimum*norm(σ_n) 371 372 σ_t = μ_d_minimum*norm(σ_n)*σ_t/norm(σ_t) 373 δ_t = σ_t*A_h/(-k_t) 374 375 E_shear = norm(σ_t*A_h)*norm(vel_lin)*simulation.time_step 376 simulation.grains[i].thermal_energy += 0.5*E_shear 377 simulation.grains[j].thermal_energy += 0.5*E_shear 378 end 379 end 380 381 force_n = dot(σ_t, n)*A_h 382 force_t = dot(σ_t, t)*A_h 383 end 384 385 # Break bond under extension through bending failure 386 if δ_n < 0.0 && tensile_strength > 0.0 && 387 shear_strength > 0.0 && norm(force_t)/A_ij > shear_strength 388 389 force_n = 0. 390 force_t = 0. 391 simulation.grains[i].contacts[ic] = 0 # remove contact 392 simulation.grains[i].n_contacts -= 1 393 simulation.grains[j].n_contacts -= 1 394 else 395 simulation.grains[i].contact_age[ic] += simulation.time_step 396 end 397 398 if simulation.grains[i].compressive_failure[ic] 399 simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t) 400 else 401 simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t.*t) 402 end 403 simulation.grains[i].contact_rotation[ic] = [0., 0., θ_t] 404 405 simulation.grains[i].force += vecTo3d(force_n.*n + force_t.*t); 406 simulation.grains[j].force -= vecTo3d(force_n.*n + force_t.*t); 407 408 simulation.grains[i].torque[3] += -force_t*R_ij + M_t 409 simulation.grains[j].torque[3] += -force_t*R_ij - M_t 410 411 simulation.grains[i].pressure += 412 force_n/simulation.grains[i].side_surface_area; 413 simulation.grains[j].pressure += 414 force_n/simulation.grains[j].side_surface_area; 415 nothing 416 end