tinteraction.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tinteraction.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