fix:remove normal force from rigid body rod contact

GazzolaLab · Jul 1, 2022 · a1c7d53 · a1c7d53
1 parent a241f65
commit a1c7d53
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 12 deletions.
diff --git a/elastica/joint.py b/elastica/joint.py
@@ -423,22 +423,20 @@ def _calculate_contact_forces_rod_rigid_body(
         # As a quick fix, use this instead
         mask = (gamma > 0.0) * 1.0
 
-        contact_force = contact_k * gamma
+        # Compute contact spring force
+        contact_force = contact_k * gamma * distance_vector
         interpenetration_velocity = (
             0.5 * (velocity_rod[..., i] + velocity_rod[..., i + 1])
             - velocity_cylinder[..., 0]
         )
-        contact_damping_force = contact_nu * _dot_product(
-            interpenetration_velocity, distance_vector
+        # Compute contact damping
+        normal_interpenetration_velocity = (
+            _dot_product(interpenetration_velocity, distance_vector) * distance_vector
         )
+        contact_damping_force = contact_nu * normal_interpenetration_velocity
 
         # magnitude* direction
-        net_contact_force = (
-            0.5
-            * mask
-            * (normal_force + (contact_damping_force + contact_force))
-            * distance_vector
-        )
+        net_contact_force = 0.5 * mask * (contact_damping_force + contact_force)
 
         # Torques acting on the cylinder
         moment_arm = x_cylinder_contact_point - x_cylinder_center

diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py
@@ -15,7 +15,7 @@ class RodCylinderParallelContact(
     rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact()
 
     # time step etc
-    final_time = 5.0
+    final_time = 10.0
     time_step = 1e-4
     total_steps = int(final_time / time_step) + 1
     rendering_fps = 30  # 20 * 1e1
@@ -28,7 +28,7 @@ class RodCylinderParallelContact(
     poisson_ratio = 0.5
     shear_modulus = E / (2 * (1 + poisson_ratio))
     n_elem = 50
-    nu = 0.0
+    nu = 0.5
     start = np.zeros((3,))
     direction = np.array([np.sin(inclination_angle), 0.0, np.cos(inclination_angle)])
     normal = np.array([0.0, 1.0, 0.0])
@@ -68,7 +68,9 @@ class RodCylinderParallelContact(
 
     # Add contact between rigid body and rod
     rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using(
-        ExternalContact, k=1e3, nu=0.1
+        ExternalContact,
+        k=5e4,
+        nu=0.1,
     )
 
     # Add callbacks