DEM improvements (#756)

* fixes

* fixing

* fix and add tests

* fix tests

* fix test

* revert change

* fix tests

* improve docs

* update news

* update

* add new example

* cleanup

* format

* Update NEWS.md

* review changes

* format

* format

* remove static arrays from dep

* review comments

* update

* remove test_util.jl

* review comments

* fix style of example

* remove again

* remove again

* update

* set 0.3

* fix pertub

Author: svchb

NEWS.md

@@ -20,6 +20,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 ### Features
 
+- **Explicit Contact Models:** Added explicit contact models, `LinearContactModel` and `HertzContactModel`, to the DEM solver. (#756)
+
 - **Particle Shifting Technique (PST) for Closed Systems:** Integrated the
   Particle Shifting Technique to enhance particle distribution, reduce clumping
   and prevent void regions due to tensile instability in closed system simulations. (#735)
@@ -88,6 +90,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - Add particle packing for 2D (.asc) and 3D (.stl) geometries (#529)
 
 ### Compatibility Changes
+
 - Dropped support for Julia 1.9
 
 ## Version 0.2.4

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
-version = "0.3-dev"
+version = "0.3"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/src/systems/dem.md

@@ -1,4 +1,5 @@
 # [Discrete Element Method](@id dem)
+
 The Discrete Element Method (DEM) is a computational technique widely used in physics, engineering,
 and applied mathematics for simulating the mechanical behavior of granular materials, such as powders,
 sand, soil, or rock, as well as other discontinua. Unlike continuum mechanics that treats materials as
@@ -7,13 +8,23 @@ detailed insights into the micro-mechanical behavior of materials, making it par
 in fields such as geomechanics, material science, and mechanical engineering.
 
 ## Fundamental Principles
+
 The core idea behind DEM is the discretization of a material system into a finite set of distinct,
 interacting mass elements (particles). These elements (particles) can vary in shape, size, and properties, and
 they interact with each other and possibly with their boundaries through contact forces and potential fields.
 The motion and behavior of each mass element are governed by Newton's laws of motion, accounting for the forces
 and moments acting upon them.
 
+## API
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "solid", "discrete_element_method", "system.jl")]
 ```
+
+### Contact Models
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "solid", "discrete_element_method", "contact_models.jl")]
+```

examples/dem/collapsing_sand_pile_3d.jl

@@ -0,0 +1,85 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+# Physical parameters
+gravity = -9.81
+acceleration = (0.0, 0.0, gravity)
+
+# ==========================================================================================
+# ==== Collapsing Sand Pile Simulation
+
+# Resolution
+particle_spacing = 0.1
+
+# Initial sand column dimensions and placement
+pile_size = (0.5, 0.5, 1.0)
+pile_min_z = particle_spacing * 0.1
+pile_center_z = pile_min_z + pile_size[3] / 2
+pile_center = (0.0, 0.0, pile_center_z)
+sand_density = 1600.0
+
+# Container dimensions
+container_width = 10
+container_depth = 10
+container_height = 1.5
+
+n_boundary_layers = 1
+boundary_thickness = n_boundary_layers * particle_spacing
+
+# Sand column particles
+n_particles_pile = round.(Int, pile_size ./ particle_spacing)
+
+min_coords_pile = (pile_center[1] - pile_size[1] / 2,
+                   pile_center[2] - pile_size[2] / 2,
+                   pile_min_z)
+sand_particles = RectangularShape(particle_spacing, n_particles_pile, min_coords_pile;
+                                  density=sand_density, coordinates_perturbation=0.1)
+
+# Boundary particles (floor and walls)
+min_boundary = (-container_width / 2, -container_depth / 2, 0.0)
+boundary_density = sand_density
+
+floor_width = container_width + 2 * boundary_thickness
+floor_depth = container_depth + 2 * boundary_thickness
+n_particles_floor_x = round(Int, floor_width / particle_spacing)
+n_particles_floor_y = round(Int, floor_depth / particle_spacing)
+n_particles_floor_z = n_boundary_layers
+
+min_coords_floor = (min_boundary[1] - boundary_thickness,
+                    min_boundary[2] - boundary_thickness,
+                    min_boundary[3] - boundary_thickness)
+floor_particles = RectangularShape(particle_spacing,
+                                   (n_particles_floor_x, n_particles_floor_y,
+                                    n_particles_floor_z),
+                                   min_coords_floor; density=boundary_density, tlsph=true)
+boundary_particles = floor_particles
+
+# ==========================================================================================
+# ==== Systems
+contact_model = LinearContactModel(1e6)
+damping_coefficient = 0.00001
+
+sand_system = DEMSystem(sand_particles, contact_model;
+                        damping_coefficient=damping_coefficient,
+                        acceleration=acceleration, radius=0.4 * particle_spacing)
+
+boundary_stiffness = 1.0e5
+boundary_system = BoundaryDEMSystem(boundary_particles, boundary_stiffness)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(sand_system, boundary_system)
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=2000)
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+callbacks = CallbackSet(info_callback, saving_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-3, # Limit stepsize to prevent crashing
+            dt=1e-7,  # Initial step size
+            save_everystep=false, callback=callbacks);

examples/dem/rectangular_tank_2d.jl

@@ -4,7 +4,12 @@ using OrdinaryDiffEq
 gravity = -9.81
 
 # ==========================================================================================
-# ==== Falling rocks
+# ==== Falling Rocks Simulation Setup
+#
+# This example sets up a simulation of falling rocks within a rectangular tank.
+# The tank is generated using a helper function (RectangularTank) which creates
+# a "fluid" region (containing the rock particles) and a "boundary" region (representing walls).
+# ==========================================================================================
 
 particle_spacing = 0.1
 
@@ -15,29 +20,48 @@ rock_density = 3000.0
 tank_width = 2.0
 tank_height = 4.0
 
+# Create a rectangular tank. The tank has a "fluid" region for the rock particles
+# and a "boundary" region for the container walls.
 tank = RectangularTank(particle_spacing, (rock_width, rock_height),
                        (tank_width, tank_height), rock_density,
                        n_layers=2)
 
 # ==========================================================================================
 # ==== Systems
+#
+# We adjust the rock positions upward to allow them to fall under gravity.
+# Next, we create a contact model and use it to build the rock (DEM) system.
+# ==========================================================================================
 
-# Move the rocks up to let them fall
+# Move the rock particles up to let them fall
 tank.fluid.coordinates[2, :] .+= 0.5
-rock_system = DEMSystem(tank.fluid, 2 * 10e5, 10e9, 0.3, acceleration=(0.0, gravity))
+# small perturbation
+tank.fluid.coordinates .+= 0.01 .* (2 .* rand(size(tank.fluid.coordinates)) .- 1)
+
+# Create a contact model.
+# Option 1: Hertzian contact model (uses elastic modulus and Poisson's ratio)
+contact_model = HertzContactModel(10e9, 0.3)
+# Option 2 (alternative): Linear contact model (constant normal stiffness)
+# contact_model = LinearContactModel(2 * 10e5)
+
+# Construct the rock system using the new DEMSystem signature.
+rock_system = DEMSystem(tank.fluid, contact_model; damping_coefficient=0.0001,
+                        acceleration=(0.0, gravity), radius=0.4 * particle_spacing)
+
+# Construct the boundary system for the tank walls.
 boundary_system = BoundaryDEMSystem(tank.boundary, 10e7)
 
 # ==========================================================================================
 # ==== Simulation
+# ==========================================================================================
 
 semi = Semidiscretization(rock_system, boundary_system)
 
-tspan = (0.0, 5.0)
+tspan = (0.0, 4.0)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=5000)
 saving_callback = SolutionSavingCallback(dt=0.02)
-
 callbacks = CallbackSet(info_callback, saving_callback)
 
 # Use a Runge-Kutta method with automatic (error based) time step size control

src/TrixiParticles.jl

@@ -79,6 +79,7 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
        PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni,
        BernoulliPressureExtrapolation
+export HertzContactModel, LinearContactModel
 export BoundaryMovement
 export examples_dir, validation_dir
 export trixi2vtk, vtk2trixi

src/io/write_vtk.jl

@@ -257,6 +257,15 @@ function write2vtk!(vtk, v, u, t, system; write_meta_data=true)
     return vtk
 end
 
+function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true)
+    vtk["velocity"] = view(v, 1:ndims(system), :)
+    vtk["mass"] = [hydrodynamic_mass(system, particle)
+                   for particle in active_particles(system)]
+    vtk["radius"] = [particle_radius(system, particle)
+                     for particle in active_particles(system)]
+    return vtk
+end
+
 function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
     vtk["velocity"] = [current_velocity(v, system, particle)
                        for particle in active_particles(system)]

src/schemes/boundary/system.jl

@@ -121,7 +121,6 @@ function Base.show(io::IO, system::BoundaryDEMSystem)
     @nospecialize system # reduce precompilation time
 
     print(io, "BoundaryDEMSystem{", ndims(system), "}(")
-    print(io, system.boundary_model)
     print(io, ") with ", nparticles(system), " particles")
 end
 
@@ -267,6 +266,10 @@ end
     return system.coordinates
 end
 
+@inline function current_velocity(v, system::BoundaryDEMSystem, particle)
+    return zero(SVector{ndims(system), eltype(system)})
+end
+
 @inline function current_velocity(v, system::BoundarySPHSystem, particle)
     return current_velocity(v, system, system.movement, particle)
 end