Granular.jl

commit 386bc4ddefe8ed8700984bc6802f5490935edbe5
parent d47ca02a8b5af731e8ed4c76d647e26a75af518c
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Fri, 26 May 2017 16:23:29 -0400

find interaction normal force in VTK writing function

Diffstat:
M
src/interaction.jl
 | 
4
++--
M
src/io.jl
 | 
24
+++++++++++++++++++++++-

2 files changed, 25 insertions(+), 3 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -76,8 +76,8 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
     δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
 
     # Effective radius
-    R_ij = harmonicMean(simulation.ice_floes[i].contact_radius,
-                        simulation.ice_floes[j].contact_radius) - abs(δ_n)/2.
+    R_ij = harmonicMean(r_i, r_j) - abs(δ_n)/2.
+
     # Contact area
     A_ij = R_ij*min(simulation.ice_floes[i].thickness, 
                     simulation.ice_floes[j].thickness)
diff --git a/src/io.jl b/src/io.jl
@@ -145,11 +145,33 @@ function writeIceFloeBondVTK(simulation::Simulation,
         for ic=1:Nc_max
             if simulation.ice_floes[i].contacts[ic] > 0
                 j = simulation.ice_floes[i].contacts[ic]
+
+                p = simulation.ice_floes[i].lin_pos -
+                    simulation.ice_floes[j].lin_pos
+                dist = norm(p)
+
+                r_i = simulation.ice_floes[i].contact_radius
+                r_j = simulation.ice_floes[j].contact_radius
+                δ_n = dist - (r_i + r_j)
+
+                if simulation.ice_floes[i].youngs_modulus > 0. &&
+                    simulation.ice_floes[j].youngs_modulus > 0.
+                    R_ij = harmonicMean(r_i, r_j)
+                    A_ij = R_ij*min(simulation.ice_floes[i].thickness, 
+                                    simulation.ice_floes[j].thickness)
+                    k_n = E*A_ij/R_ij
+                else
+                    k_n = harmonicMean(simulation.ice_floes[i].
+                                       contact_stiffness_normal,
+                                       simulation.ice_floes[j].
+                                       contact_stiffness_normal)
+                end
+
                 
                 append!(i1, i)
                 append!(i2, j)
 
-                append!(force, f_n)
+                append!(force, k_n*δ_n)
 
                 append!(shear_displacement_1, simulation.ice_floes[i].
                         contact_parallel_displacement[ic][1])