Rigid contact friction history

docs/literate/src/tut_rigid_body_fsi.jl

@@ -227,10 +227,15 @@ nothing # hide
 
 # ## Step 3: With contact model
 # Finally, we add a `contact_model` to handle collisions between rigid bodies and between
-# rigid bodies and the tank.
+# rigid bodies and the tank. Here we also enable the frictional rigid-wall path, so we need
+# `UpdateCallback()` to update the tangential contact history between time steps.
 
 contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                   normal_damping=150.0,
+                                  static_friction_coefficient=0.6,
+                                  kinetic_friction_coefficient=0.4,
+                                  tangential_stiffness=1.0e5,
+                                  tangential_damping=150.0,
                                   contact_distance=2.0 * structure_particle_spacing)
 nothing # hide
 
@@ -253,15 +258,15 @@ nothing # hide
 
 info_callback = InfoCallback(interval=100) # hide
 saving_callback_step3 = SolutionSavingCallback(dt=0.02, prefix="step3") # hide
-callbacks_step3 = CallbackSet(info_callback, saving_callback_step3) # hide
+callbacks_step3 = CallbackSet(info_callback, saving_callback_step3, UpdateCallback()) # hide
 nothing # hide
 
 # ```julia
 # sol_step3 = solve(ode_step3, RDPK3SpFSAL49(), save_everystep=false,
-#                   callback=callbacks,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
+#                   callback=callbacks_step3,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
 # ```
 sol_step3 = solve(ode_step3, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                  save_everystep=false) # hide
+                  save_everystep=false, callback=callbacks_step3) # hide
 nothing # hide
 
 # The plot now shows the full simulation. The squares collide with the tank bottom and each other.
@@ -320,15 +325,15 @@ nothing # hide
 
 info_callback = InfoCallback(interval=100) # hide
 saving_callback_step4 = SolutionSavingCallback(dt=0.02, prefix="step4") # hide
-callbacks_step4 = CallbackSet(info_callback, saving_callback_step4) # hide
+callbacks_step4 = CallbackSet(info_callback, saving_callback_step4, UpdateCallback()) # hide
 nothing # hide
 
 # ```julia
 # sol_step4 = solve(ode_step4, RDPK3SpFSAL49(), save_everystep=false,
-#                   callback=callbacks,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
+#                   callback=callbacks_step4,  abstol=1e-6, reltol=1e-4, dtmax=2e-3)
 # ```
 sol_step4 = solve(ode_step4, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                  save_everystep=false) # hide
+                  save_everystep=false, callback=callbacks_step4) # hide
 nothing # hide
 
 # And here is the final plot with circles instead of squares.
@@ -372,9 +377,10 @@ semi_next = Semidiscretization(fluid_system, boundary_system,
                                small_sphere_systems...)
 
 ode_step_next = semidiscretize(semi_next, tspan) #hide
+callbacks_step_next = CallbackSet(UpdateCallback()) # hide
 sol_step_next = solve(ode_step_next, RDPK3SpFSAL49(), # hide
                       abstol=1e-6, reltol=1e-4, dtmax=2e-3, # hide
-                      save_everystep=false) # hide
+                      save_everystep=false, callback=callbacks_step_next) # hide
 nothing # hide
 
 plot(sol_step_next, legend=nothing) #hide
@@ -428,8 +434,9 @@ hexagon_system = RigidBodySystem(hexagon_shape;
 semi_hexagon = Semidiscretization(fluid_system, boundary_system, hexagon_system)
 ode_step_hex = semidiscretize(semi_hexagon, (0.0, 0.4))
 
+callbacks_step_hex = CallbackSet(UpdateCallback()) # hide
 sol_step_hex = solve(ode_step_hex, RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-5, dtmax=1e-3,
-                     save_everystep=false)
+                     save_everystep=false, callback=callbacks_step_hex)
 plot(sol_step_hex, legend=nothing) #hide
 plot!(dpi=200) # hide
 savefig("tut_rigid_body_fsi_hex.png"); # hide

docs/src/callbacks.md

@@ -1,5 +1,10 @@
 # Callbacks
 
+[`UpdateCallback`](@ref) is required for systems that keep mutable state between time
+steps. In the current rigid-contact implementation, this applies when a
+[`RigidContactModel`](@ref) uses tangential history, i.e. the rigid-wall friction path
+needs the tangential displacement cache to be updated between steps.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = map(file -> joinpath("callbacks", file), readdir(joinpath("..", "src", "callbacks")))

docs/src/systems/rigid_body.md

@@ -17,6 +17,39 @@ Rigid contact is configured through the contact model. This is separate from the
 boundary model used for fluid-structure interaction; see
 [Boundary Models](@ref boundary_models) for that part of the rigid-body setup.
 
+`RigidContactModel` defines the rigid-contact law shared by rigid-wall and rigid-rigid
+interaction. The always-active parameters are `normal_stiffness`, `normal_damping`, and
+`contact_distance`.
+
+For rigid-wall contact, the current implementation also supports tangential friction with
+the parameters `static_friction_coefficient`, `kinetic_friction_coefficient`,
+`tangential_stiffness`, `tangential_damping`, `stick_velocity_tolerance`, and
+`penetration_slop`. When the tangential spring history is active, this requires
+[`UpdateCallback`](@ref) so the tangential displacement cache is updated between time
+steps.
+
+The current implementation uses the same model for rigid-wall and rigid-rigid contact, but
+the active behavior is not identical yet:
+
+- rigid-wall contact groups penetrating wall neighbors into a small number of contact
+  manifolds per rigid particle and applies one normal-plus-tangential contact force per
+  manifold,
+- rigid-rigid contact evaluates direct pairwise normal contact forces between rigid
+  particles,
+- and rigid-rigid remains normal-only for now, even though the tangential-history key
+  design is already shared with the wall path.
+
+`contact_distance` defines when contact starts. If `contact_distance == 0`, the
+particle spacing of the `RigidBodySystem` is used when the contact model is adapted to
+the runtime system.
+
+If no `contact_model` is specified for a rigid body, rigid-wall and rigid-rigid contact
+for that system are disabled.
+
+For output and postprocessing, rigid bodies also expose the diagnostics
+`contact_count` and `max_contact_penetration`. They are available through rigid-body
+system data and VTK output.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "structure", "rigid_body", "contact_models.jl")]

examples/fsi/falling_rigid_spheres_2d.jl

@@ -99,9 +99,14 @@ boundary_model_structure_2 = BoundaryModelDummyParticles(hydrodynamic_densities_
                                                          fluid_smoothing_kernel,
                                                          fluid_smoothing_length)
 
-# Basic rigid contact model used for both rigid bodies.
+# Use the frictional rigid-wall contact path, which requires `UpdateCallback()` to keep
+# the tangential contact history in sync between time steps.
 contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                   normal_damping=200.0,
+                                  static_friction_coefficient=0.6,
+                                  kinetic_friction_coefficient=0.4,
+                                  tangential_stiffness=1.0e5,
+                                  tangential_damping=150.0,
                                   contact_distance=2.0 *
                                                    structure_particle_spacing)
 
@@ -125,7 +130,7 @@ ode = semidiscretize(semi, tspan)
 info_callback = InfoCallback(interval=50)
 saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="")
 
-callbacks = CallbackSet(info_callback, saving_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 sol = solve(ode, RDPK3SpFSAL49(),

examples/fsi/falling_rotating_rigid_squares_2d.jl

@@ -109,13 +109,22 @@ end
 boundary_model_structure_1 = structure_boundary_model(square1)
 boundary_model_structure_2 = structure_boundary_model(square2)
 
-# Use a less dissipative wall contact for the denser square so its rebound is more visible.
+# Use frictional rigid-wall contact so the rotating squares exchange tangential impulses with
+# the tank floor as well as normal contact forces. This requires `UpdateCallback()`.
 contact_model_1 = RigidContactModel(; normal_stiffness=2.0e5,
                                     normal_damping=200.0,
+                                    static_friction_coefficient=0.6,
+                                    kinetic_friction_coefficient=0.4,
+                                    tangential_stiffness=1.0e5,
+                                    tangential_damping=180.0,
                                     contact_distance=2.0 *
                                                      structure_particle_spacing)
 contact_model_2 = RigidContactModel(; normal_stiffness=2.0e5,
                                     normal_damping=80.0,
+                                    static_friction_coefficient=0.5,
+                                    kinetic_friction_coefficient=0.3,
+                                    tangential_stiffness=8.0e4,
+                                    tangential_damping=120.0,
                                     contact_distance=2.0 *
                                                      structure_particle_spacing)
 
@@ -144,7 +153,7 @@ saving_callback = SolutionSavingCallback(dt=0.01,
                                          output_directory="out",
                                          prefix="")
 
-callbacks = CallbackSet(info_callback, saving_callback)
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 # To prevent penetration of fluid particles through the rigid bodies or the boundary

examples/fsi/falling_rotating_rigid_squares_w_buoys_2d.jl

@@ -20,6 +20,10 @@ small_sphere_y = initial_fluid_size[2] + small_sphere_radius
 small_sphere_x_positions = 0.2:(3 * small_sphere_radius):1.8
 small_sphere_contact_model = RigidContactModel(; normal_stiffness=2.0e5,
                                                normal_damping=120.0,
+                                               static_friction_coefficient=0.6,
+                                               kinetic_friction_coefficient=0.4,
+                                               tangential_stiffness=1.0e5,
+                                               tangential_damping=150.0,
                                                contact_distance=2.0 *
                                                                 structure_particle_spacing)
 extra_structure_systems = [begin

examples/structure/sliding_rigid_squares_friction_2d.jl

@@ -0,0 +1,99 @@
+# ==========================================================================================
+# 2D Sliding Rigid Squares with and without Wall Friction
+#
+# Two identical rigid squares slide on the same floor. The left square uses the normal-only
+# rigid contact model from PR1, while the right square uses the frictional wall-contact path
+# added in PR2. The frictional square slows down and starts rotating due to tangential wall
+# forces, whereas the normal-only square keeps sliding without spin-up.
+#
+# In ParaView, compare the trajectories and the rigid-body field data such as
+# `angular_velocity`, `contact_count`, and `max_contact_penetration`.
+# ==========================================================================================
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.03
+boundary_layers = 3
+
+# ==========================================================================================
+# ==== Experiment Setup
+gravity = 9.81
+tspan = (0.0, 0.8)
+
+square_side_length = 0.18
+square_density = 1000.0
+square_particles_per_side = round(Int, square_side_length / particle_spacing)
+square_bottom_y = 0.03
+
+square_frictionless = RectangularShape(particle_spacing,
+                                       (square_particles_per_side,
+                                        square_particles_per_side),
+                                       (-1.0, square_bottom_y),
+                                       density=square_density,
+                                       velocity=(1.0, 0.0))
+square_frictional = RectangularShape(particle_spacing,
+                                     (square_particles_per_side,
+                                      square_particles_per_side),
+                                     (0.55, square_bottom_y),
+                                     density=square_density,
+                                     velocity=(1.0, 0.0))
+
+# ==========================================================================================
+# ==== Wall Boundary
+floor_length = 3.0
+floor_height = 0.03
+wall_density = 1000.0
+
+floor = RectangularTank(particle_spacing, (0.0, 0.0), (floor_length, floor_height),
+                        wall_density, n_layers=boundary_layers,
+                        min_coordinates=(-1.5, 0.0),
+                        faces=(false, false, true, false))
+
+boundary_model = BoundaryModelMonaghanKajtar(10.0, 1.0, particle_spacing,
+                                             floor.boundary.mass)
+boundary_system = WallBoundarySystem(floor.boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Rigid Structures
+contact_model_frictionless = RigidContactModel(; normal_stiffness=2.0e5,
+                                               normal_damping=180.0,
+                                               contact_distance=2.0 * particle_spacing)
+
+contact_model_frictional = RigidContactModel(; normal_stiffness=2.0e5,
+                                             normal_damping=180.0,
+                                             static_friction_coefficient=0.6,
+                                             kinetic_friction_coefficient=0.4,
+                                             tangential_stiffness=1.0e5,
+                                             tangential_damping=150.0,
+                                             contact_distance=2.0 * particle_spacing)
+
+structure_system_frictionless = RigidBodySystem(square_frictionless;
+                                                contact_model=contact_model_frictionless,
+                                                acceleration=(0.0, -gravity),
+                                                particle_spacing=particle_spacing,
+                                                color_value=1)
+structure_system_frictional = RigidBodySystem(square_frictional;
+                                              contact_model=contact_model_frictional,
+                                              acceleration=(0.0, -gravity),
+                                              particle_spacing=particle_spacing,
+                                              color_value=2)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(structure_system_frictionless, structure_system_frictional,
+                          boundary_system)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-6,
+            reltol=1e-4,
+            dtmax=1e-3,
+            save_everystep=false, callback=callbacks);

src/callbacks/post_process.jl

@@ -159,7 +159,8 @@ function (pp::PostprocessCallback)(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
+                                                     reset_interaction_caches=false)
         end
 
         # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.

src/callbacks/update.jl

@@ -96,6 +96,10 @@ function (update_callback!::UpdateCallback)(integrator)
             update_particle_packing(system, v_ode, u_ode, semi, integrator)
         end
 
+        foreach_system(semi) do system
+            update_rigid_contact_eachstep!(system, v_ode, u_ode, semi, t, integrator)
+        end
+
         # This is only used by the particle packing system and should be removed in the future
         foreach_system(semi) do system
             update_transport_velocity!(system, v_ode, semi, integrator)

src/general/abstract_system.jl

@@ -153,7 +153,7 @@ function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_final!(system, v, u, v_ode, u_ode, semi, t)
+function update_final!(system, v, u, v_ode, u_ode, semi, t; kwargs...)
     return system
 end
 

src/general/semidiscretization.jl

@@ -491,7 +491,7 @@ end
 
 # Update the systems and neighborhood searches (NHS) for a simulation
 # before calling `interact!` to compute forces.
-function update_systems_and_nhs(v_ode, u_ode, semi, t)
+function update_systems_and_nhs(v_ode, u_ode, semi, t; reset_interaction_caches=true)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the TLSPH system)
     foreach_system(semi) do system
@@ -539,7 +539,8 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
-        update_final!(system, v, u, v_ode, u_ode, semi, t)
+        update_final!(system, v, u, v_ode, u_ode, semi, t;
+                      reset_interaction_caches)
     end
 end
 

src/io/io.jl

@@ -128,6 +128,7 @@ function add_system_data!(system_data, system::RigidBodySystem)
     system_data["adhesion_coefficient"] = system.adhesion_coefficient
     system_data["color"] = system.cache.color
     add_system_data!(system_data, system.boundary_model)
+    add_system_data!(system_data, system.contact_model)
 end
 
 function add_system_data!(system_data, system::WallBoundarySystem)
@@ -219,6 +220,20 @@ function add_system_data!(system_data, contact_model::LinearContactModel)
     system_data["contact_model"]["normal_stiffness"] = contact_model.normal_stiffness
 end
 
+function add_system_data!(system_data, contact_model::RigidContactModel)
+    system_data["contact_model"] = Dict{String, Any}()
+    system_data["contact_model"]["model"] = type2string(contact_model)
+    system_data["contact_model"]["normal_stiffness"] = contact_model.normal_stiffness
+    system_data["contact_model"]["normal_damping"] = contact_model.normal_damping
+    system_data["contact_model"]["static_friction_coefficient"] = contact_model.static_friction_coefficient
+    system_data["contact_model"]["kinetic_friction_coefficient"] = contact_model.kinetic_friction_coefficient
+    system_data["contact_model"]["tangential_stiffness"] = contact_model.tangential_stiffness
+    system_data["contact_model"]["tangential_damping"] = contact_model.tangential_damping
+    system_data["contact_model"]["contact_distance"] = contact_model.contact_distance
+    system_data["contact_model"]["stick_velocity_tolerance"] = contact_model.stick_velocity_tolerance
+    system_data["contact_model"]["penetration_slop"] = contact_model.penetration_slop
+end
+
 function add_system_data!(system_data, state_equation::StateEquationCole)
     system_data["state_equation"] = Dict{String, Any}()
     system_data["state_equation"]["model"] = type2string(state_equation)