commit 8268169fbd4edcf11db3621ece0e1496d9b0e25d
parent 5ac16e5a720ebb6d1d93d253f3d8a06ab2800b3a
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Tue, 13 Jul 2021 14:58:58 +0200
include tangential interactions with non-rotating particles
to get the old rotate=false behavior physics, set both:
grain.contact_dynamic_friction = 0.0
grain.contact_viscosity_tangential = 0.0
Diffstat:
| M | src/interaction.jl | | | 120 | ++++++++++++++++++++++++++++++------------------------------------------------- |
1 file changed, 45 insertions(+), 75 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -145,31 +145,23 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.grains[j].lin_vel[1:2]
vel_n = dot(vel_lin, n)
- if !simulation.grains[i].rotating && !simulation.grains[j].rotating
- rotation = false
+ vel_t = dot(vel_lin, t) -
+ harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel[3] +
+ simulation.grains[j].ang_vel[3])
+
+ # Correct old tangential displacement for contact rotation and add new
+ δ_t0 = simulation.grains[i].contact_parallel_displacement[ic][1:2]
+ if !simulation.grains[i].compressive_failure[ic]
+ δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
else
- rotation = true
+ # Use δ_t as total slip-distance vector
+ δ_t = δ_t0 .+ (vel_lin .+ vel_t.*t) .* simulation.time_step
end
-
- if rotation
- vel_t = dot(vel_lin, t) -
- harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel[3] +
- simulation.grains[j].ang_vel[3])
-
- # Correct old tangential displacement for contact rotation and add new
- δ_t0 = simulation.grains[i].contact_parallel_displacement[ic][1:2]
- if !simulation.grains[i].compressive_failure[ic]
- δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
- else
- # Use δ_t as total slip-distance vector
- δ_t = δ_t0 .+ (vel_lin .+ vel_t.*t) .* simulation.time_step
- end
- # Determine the contact rotation (2d)
- θ_t = simulation.grains[i].contact_rotation[ic][3] +
- (simulation.grains[j].ang_vel[3] -
- simulation.grains[i].ang_vel[3]) * simulation.time_step
- end
+ # Determine the contact rotation (2d)
+ θ_t = simulation.grains[i].contact_rotation[ic][3] +
+ (simulation.grains[j].ang_vel[3] -
+ simulation.grains[i].ang_vel[3]) * simulation.time_step
# Effective radius
R_ij = harmonicMean(r_i, r_j)
@@ -198,30 +190,24 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Effective normal and tangential stiffness
k_n = E * A_ij/R_ij
#k_t = k_n*ν # Kneib et al 2016
- if rotation
- k_t = k_n * 2. * (1. - ν^2.) / ((2. - ν) * (1. + ν)) # Obermayr 2011
- end
+ k_t = k_n * 2. * (1. - ν^2.) / ((2. - ν) * (1. + ν)) # Obermayr 2011
else # Micromechanical parameterization: k_n and k_t set explicitly
k_n = harmonicMean(simulation.grains[i].contact_stiffness_normal,
simulation.grains[j].contact_stiffness_normal)
- if rotation
- k_t = harmonicMean(simulation.grains[i].contact_stiffness_tangential,
- simulation.grains[j].contact_stiffness_tangential)
- end
+ k_t = harmonicMean(simulation.grains[i].contact_stiffness_tangential,
+ simulation.grains[j].contact_stiffness_tangential)
end
γ_n = harmonicMean(simulation.grains[i].contact_viscosity_normal,
simulation.grains[j].contact_viscosity_normal)
- if rotation
- γ_t = harmonicMean(simulation.grains[i].contact_viscosity_tangential,
- simulation.grains[j].contact_viscosity_tangential)
+ γ_t = harmonicMean(simulation.grains[i].contact_viscosity_tangential,
+ simulation.grains[j].contact_viscosity_tangential)
- μ_d_minimum = min(simulation.grains[i].contact_dynamic_friction,
- simulation.grains[j].contact_dynamic_friction)
- end
+ μ_d_minimum = min(simulation.grains[i].contact_dynamic_friction,
+ simulation.grains[j].contact_dynamic_friction)
# Determine contact forces
if k_n ≈ 0. && γ_n ≈ 0. # No interaction
@@ -257,39 +243,28 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
tensile_strength = min(simulation.grains[i].contact_age[ic]*
simulation.grains[i].strength_heal_rate,
simulation.grains[i].tensile_strength)
- if rotation
shear_strength = min(simulation.grains[i].contact_age[ic]*
simulation.grains[i].strength_heal_rate,
simulation.grains[i].shear_strength)
- M_t = 0.0
- if tensile_strength > 0.0
- # Determine bending momentum on contact [N*m],
- # (converting k_n to E to bar(k_n))
- M_t = (k_n*R_ij/(A_ij*(simulation.grains[i].contact_radius +
+ M_t = 0.0
+ if tensile_strength > 0.0
+ # Determine bending momentum on contact [N*m],
+ # (converting k_n to E to bar(k_n))
+ M_t = (k_n*R_ij/(A_ij*(simulation.grains[i].contact_radius +
simulation.grains[j].contact_radius)))*I_ij*θ_t
- end
end
# Reset contact age (breaking bond) if bond strength is exceeded
- if rotation
- if δ_n >= 0.0 && norm(force_n)/A_ij + norm(M_t)*R_ij/I_ij > tensile_strength
- force_n = 0.
- force_t = 0.
- simulation.grains[i].contacts[ic] = 0 # remove contact
- simulation.grains[i].n_contacts -= 1
- simulation.grains[j].n_contacts -= 1
- end
- else
- if δ_n >= 0.0 && norm(force_n)/A_ij > tensile_strength
- force_n = 0.
- force_t = 0.
- simulation.grains[i].contacts[ic] = 0 # remove contact
- simulation.grains[i].n_contacts -= 1
- simulation.grains[j].n_contacts -= 1
- end
+
+ if δ_n >= 0.0 && norm(force_n)/A_ij + norm(M_t)*R_ij/I_ij > tensile_strength
+ force_n = 0.
+ force_t = 0.
+ simulation.grains[i].contacts[ic] = 0 # remove contact
+ simulation.grains[i].n_contacts -= 1
+ simulation.grains[j].n_contacts -= 1
end
- if rotation && !simulation.grains[i].compressive_failure[ic]
+ if !simulation.grains[i].compressive_failure[ic]
if k_t ≈ 0. && γ_t ≈ 0.
# do nothing
@@ -408,7 +383,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
# Break bond under extension through bending failure
- if δ_n < 0.0 && tensile_strength > 0.0 && rotation &&
+ if δ_n < 0.0 && tensile_strength > 0.0 &&
shear_strength > 0.0 && norm(force_t)/A_ij > shear_strength
force_n = 0.
@@ -420,23 +395,18 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.grains[i].contact_age[ic] += simulation.time_step
end
- if rotation
- if simulation.grains[i].compressive_failure[ic]
- simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t)
- else
- simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t.*t)
- end
- simulation.grains[i].contact_rotation[ic] = [0., 0., θ_t]
-
- simulation.grains[i].force += vecTo3d(force_n.*n + force_t.*t);
- simulation.grains[j].force -= vecTo3d(force_n.*n + force_t.*t);
-
- simulation.grains[i].torque[3] += -force_t*R_ij + M_t
- simulation.grains[j].torque[3] += -force_t*R_ij - M_t
+ if simulation.grains[i].compressive_failure[ic]
+ simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t)
else
- simulation.grains[i].force += vecTo3d(force_n.*n);
- simulation.grains[j].force -= vecTo3d(force_n.*n);
+ simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t.*t)
end
+ simulation.grains[i].contact_rotation[ic] = [0., 0., θ_t]
+
+ simulation.grains[i].force += vecTo3d(force_n.*n + force_t.*t);
+ simulation.grains[j].force -= vecTo3d(force_n.*n + force_t.*t);
+
+ simulation.grains[i].torque[3] += -force_t*R_ij + M_t
+ simulation.grains[j].torque[3] += -force_t*R_ij - M_t
simulation.grains[i].pressure +=
force_n/simulation.grains[i].side_surface_area;
