From 0283d166921483126d1b5a5ac21f992433208499 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 24 Oct 2023 11:55:00 +0200 Subject: [PATCH 1/4] implement basic dem --- examples/dem/rectangular_tank_2d.jl | 66 ++++++++++++++++++ src/TrixiParticles.jl | 2 +- src/general/initial_condition.jl | 5 +- src/general/semidiscretization.jl | 6 ++ src/schemes/schemes.jl | 2 + .../discrete_element_method.jl | 1 + .../solid/discrete_element_method/rhs.jl | 38 +++++++++++ .../solid/discrete_element_method/system.jl | 68 +++++++++++++++++++ src/visualization/write2vtk.jl | 11 ++- 9 files changed, 196 insertions(+), 3 deletions(-) create mode 100644 examples/dem/rectangular_tank_2d.jl create mode 100644 src/schemes/solid/discrete_element_method/discrete_element_method.jl create mode 100644 src/schemes/solid/discrete_element_method/rhs.jl create mode 100644 src/schemes/solid/discrete_element_method/system.jl diff --git a/examples/dem/rectangular_tank_2d.jl b/examples/dem/rectangular_tank_2d.jl new file mode 100644 index 000000000..d23f7fc78 --- /dev/null +++ b/examples/dem/rectangular_tank_2d.jl @@ -0,0 +1,66 @@ +using TrixiParticles +using OrdinaryDiffEq + +gravity = -9.81 + +# ========================================================================================== +# ==== Fluid + +particle_spacing = 0.1 + +# Ratio of fluid particle spacing to boundary particle spacing +beta = 1 +boundary_layers = 1 + +water_width = 2.0 +water_height = 2.0 +water_density = 1000.0 + +tank_width = 2.0 +tank_height = 4.0 + +tank = RectangularTank(particle_spacing, (water_width, water_height), + (tank_width, tank_height), water_density, + n_layers=boundary_layers, spacing_ratio=beta) + +# ========================================================================================== +# ==== Boundary models + +# boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +# state_equation=state_equation, +# AdamiPressureExtrapolation(), +# smoothing_kernel, smoothing_length) + +# ========================================================================================== +# ==== Systems + +solid_system = DEMSystem(tank.fluid, 0.1, acceleration=(0.0, gravity)) + +# boundary_system = BoundaryDEMSystem(tank.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation + +semi = Semidiscretization(solid_system, neighborhood_search=GridNeighborhoodSearch) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +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. +# Enable threading of the RK method for better performance on multiple threads. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces +# become extremely large when fluid particles are very close to boundary particles, +# and the time integration method interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL49(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 1e1ef1508..8737c1e10 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -34,7 +34,7 @@ include("visualization/write2vtk.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem + BoundarySPHSystem, DEMSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 306a08aa1..893a16623 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -5,6 +5,7 @@ struct InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} density :: Array{ELTYPE, 1} pressure :: Array{ELTYPE, 1} + radius :: Array{ELTYPE, 1} function InitialCondition(coordinates, velocities, masses, densities; pressure=0.0, particle_spacing=-1.0) @@ -24,8 +25,10 @@ struct InitialCondition{ELTYPE} "same length as `masses`")) end + radius = 0.5 * particle_spacing * ones(length(masses)) + return new{eltype(coordinates)}(particle_spacing, coordinates, velocities, masses, - densities, pressure) + densities, pressure, radius) end end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index d91d0f8b2..48006a3f4 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -420,6 +420,12 @@ function system_interaction!(dv_ode, v_ode, u_ode, semi) end # NHS updates + +# Default to update +function nhs_coords(system, neighbor, u) + return current_coordinates(u, neighbor) +end + function nhs_coords(system::FluidSystem, neighbor::FluidSystem, u) return current_coordinates(u, neighbor) diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index 68c10c017..02e169c10 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -3,9 +3,11 @@ include("fluid/fluid.jl") include("boundary/boundary.jl") include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") +include("solid/discrete_element_method/discrete_element_method.jl") # Include rhs for all schemes include("fluid/weakly_compressible_sph/rhs.jl") include("fluid/entropically_damped_sph/rhs.jl") include("boundary/rhs.jl") include("solid/total_lagrangian_sph/rhs.jl") +include("solid/discrete_element_method/rhs.jl") diff --git a/src/schemes/solid/discrete_element_method/discrete_element_method.jl b/src/schemes/solid/discrete_element_method/discrete_element_method.jl new file mode 100644 index 000000000..c67930ab8 --- /dev/null +++ b/src/schemes/solid/discrete_element_method/discrete_element_method.jl @@ -0,0 +1 @@ +include("system.jl") diff --git a/src/schemes/solid/discrete_element_method/rhs.jl b/src/schemes/solid/discrete_element_method/rhs.jl new file mode 100644 index 000000000..92440c9ab --- /dev/null +++ b/src/schemes/solid/discrete_element_method/rhs.jl @@ -0,0 +1,38 @@ +# Calculate the interaction forces between particles in a Discrete Element Method (DEM) system. +# +# This function loops over all pairs of particles and their neighbors within a set distance. +# When particles overlap (i.e., they come into contact), a normal force is applied to resolve the overlap. +# This force is computed based on Hertzian contact mechanics typical for DEM simulations. +# The force is proportional to the amount of overlap and is directed along the normal between the particle centers. +# The magnitude of the force is determined by the stiffness constant `kn` and the overlap distan +function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, + u_neighbor_system, neighborhood_search, particle_system::DEMSystem, + neighbor_system::DEMSystem) + (; mass, radius, kn) = particle_system + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Calculate the overlap (penetration depth) between the two particles + overlap = radius[particle] + radius[neighbor] - distance + # If there's no overlap, no force needs to be applied + overlap <= 0 && return + + # Normal direction from particle to neighbor + normal = pos_diff / distance + + # Compute the force magnitude using Hertzian contact mechanics + force_magnitude = kn * overlap + + force = force_magnitude * normal + + # Update the acceleration of the particle based on the force and its mass + @inbounds for i in 1:ndims(particle_system) + dv[i, particle] += force[i] / mass[particle] + end + end + + return dv +end diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl new file mode 100644 index 000000000..cc666cf3d --- /dev/null +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -0,0 +1,68 @@ +struct DEMSystem{NDIMS, ELTYPE <: Real} <: System{NDIMS} + initial_condition :: InitialCondition{ELTYPE} + mass :: Array{ELTYPE, 1} # [particle] + radius :: Array{ELTYPE, 1} # [particle] + kn :: ELTYPE # Normal stiffness + acceleration :: SVector{NDIMS, ELTYPE} + + function DEMSystem(initial_condition, kn; acceleration=ntuple(_ -> 0.0, + ndims(initial_condition))) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + mass = copy(initial_condition.mass) + radius = copy(initial_condition.radius) + + # Make acceleration an SVector + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end + + + return new{NDIMS, ELTYPE}(initial_condition, mass, radius, kn, acceleration_) + end +end + +function Base.show(io::IO, system::DEMSystem) + @nospecialize system # reduce precompilation time + + print(io, "DEMSystem{", ndims(system), "}(") + print(io, system.initial_condition) + print(io, ", ", system.kn) + print(io, ") with ", TrixiParticles.nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::DEMSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + TrixiParticles.summary_header(io, "DEMSystem{$(ndims(system))}") + TrixiParticles.summary_line(io, "#particles", TrixiParticles.nparticles(system)) + TrixiParticles.summary_line(io, "kn", system.kn) + TrixiParticles.summary_footer(io) + end +end + +timer_name(::DEMSystem) = "solid" + +function TrixiParticles.write_u0!(u0, system::DEMSystem) + u0 .= system.initial_condition.coordinates + return u0 +end + +function TrixiParticles.write_v0!(v0, system::DEMSystem) + v0 .= system.initial_condition.velocity + return v0 +end + +# Nothing to initialize for this system +initialize!(system::DEMSystem, neighborhood_search) = system + +function compact_support(system::DEMSystem, neighbor::DEMSystem) + # we for now assume that the compact support is 3 * radius + # todo: needs to be changed for more complex simulations + return 3 * max(maximum(system.radius), maximum(neighbor.radius)) +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 8bedd5642..e9dcb9302 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -68,6 +68,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] vtk_grid(file, points, cells) do vtk + # dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) # Store particle index @@ -135,10 +136,18 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo return file end +# default to "particle" +vtkname(system) = "particle" vtkname(system::FluidSystem) = "fluid" -vtkname(system::TotalLagrangianSPHSystem) = "solid" +vtkname(system::Union{TotalLagrangianSPHSystem, DEMSystem}) = "solid" vtkname(system::BoundarySPHSystem) = "boundary" +function write2vtk!(vtk, v, u, t, system; write_meta_data=true) + vtk["velocity"] = view(v, 1:ndims(system), :) + + return vtk +end + function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["velocity"] = view(v, 1:ndims(system), :) vtk["density"] = [particle_density(v, system, particle) From b8c81405d6364b845bfa83a36d6f40d44f9bba55 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 24 Oct 2023 17:38:16 +0200 Subject: [PATCH 2/4] fix --- examples/dem/rectangular_tank_2d.jl | 2 +- src/schemes/solid/discrete_element_method/rhs.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/dem/rectangular_tank_2d.jl b/examples/dem/rectangular_tank_2d.jl index d23f7fc78..4e1d8334d 100644 --- a/examples/dem/rectangular_tank_2d.jl +++ b/examples/dem/rectangular_tank_2d.jl @@ -6,7 +6,7 @@ gravity = -9.81 # ========================================================================================== # ==== Fluid -particle_spacing = 0.1 +particle_spacing = 0.25 # Ratio of fluid particle spacing to boundary particle spacing beta = 1 diff --git a/src/schemes/solid/discrete_element_method/rhs.jl b/src/schemes/solid/discrete_element_method/rhs.jl index 92440c9ab..b6464c625 100644 --- a/src/schemes/solid/discrete_element_method/rhs.jl +++ b/src/schemes/solid/discrete_element_method/rhs.jl @@ -15,8 +15,12 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + # Calculate the overlap (penetration depth) between the two particles overlap = radius[particle] + radius[neighbor] - distance + # If there's no overlap, no force needs to be applied overlap <= 0 && return From 99b9cfdfbe3a8acead9a32eb226a6bdcf1168e69 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 25 Oct 2023 11:15:49 +0200 Subject: [PATCH 3/4] fix bnd --- examples/dem/rectangular_tank_2d.jl | 5 +--- examples/fluid/rectangular_tank_2d.jl | 5 ++-- .../dummy_particles/dummy_particles.jl | 28 +++++++++++++++---- 3 files changed, 27 insertions(+), 11 deletions(-) diff --git a/examples/dem/rectangular_tank_2d.jl b/examples/dem/rectangular_tank_2d.jl index 4e1d8334d..b9df41024 100644 --- a/examples/dem/rectangular_tank_2d.jl +++ b/examples/dem/rectangular_tank_2d.jl @@ -26,10 +26,7 @@ tank = RectangularTank(particle_spacing, (water_width, water_height), # ========================================================================================== # ==== Boundary models -# boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, -# state_equation=state_equation, -# AdamiPressureExtrapolation(), -# smoothing_kernel, smoothing_length) +boundary_model = BoundaryModelDummyParticles(tank.boundary.radius) # ========================================================================================== # ==== Systems diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl index c8fb0329c..55f318ba9 100644 --- a/examples/fluid/rectangular_tank_2d.jl +++ b/examples/fluid/rectangular_tank_2d.jl @@ -24,7 +24,8 @@ state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0, background_pressure=100000.0, clip_negative_pressure=false) -smoothing_length = 1.2 * particle_spacing +smoothing_factor = 1.2 +smoothing_length = smoothing_factor * particle_spacing smoothing_kernel = SchoenbergCubicSplineKernel{2}() viscosity = ArtificialViscosityMonaghan(0.02, 0.0) @@ -39,7 +40,7 @@ tank = RectangularTank(particle_spacing, (water_width, water_height), boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length) + smoothing_kernel, particle_spacing, smoothing_factor) # K = 9.81 * water_height # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta, diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 747b55b45..9d0aba397 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -68,28 +68,46 @@ We provide three options to compute the boundary density and pressure, determine struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C} pressure :: Vector{ELTYPE} hydrodynamic_mass :: Vector{ELTYPE} + radius :: Vector{ELTYPE} state_equation :: SE density_calculator :: DC smoothing_kernel :: K smoothing_length :: ELTYPE + smoothing_factor :: ELTYPE viscosity :: V cache :: C function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, - density_calculator, smoothing_kernel, - smoothing_length; viscosity=NoViscosity(), + density_calculator, smoothing_kernel, particle_spacing, + smoothing_factor; viscosity=NoViscosity(), state_equation=nothing) + ELTYPE = eltype(initial_density) pressure = similar(initial_density) n_particles = length(initial_density) + radius = 0.5 * ones(ELTYPE, n_particles) * particle_spacing + smoothing_length = smoothing_factor * particle_spacing + cache = (; create_cache(viscosity, n_particles, ndims(smoothing_kernel))..., create_cache(initial_density, density_calculator)...) - new{eltype(initial_density), typeof(state_equation), + new{ELTYPE, typeof(state_equation), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(cache)}(pressure, hydrodynamic_mass, state_equation, density_calculator, - smoothing_kernel, smoothing_length, viscosity, cache) + typeof(cache)}(pressure, hydrodynamic_mass, radius, state_equation, density_calculator, + smoothing_kernel, smoothing_length, smoothing_factor, viscosity, cache) + end + + function BoundaryModelDummyParticles(radius) + ELTYPE = eltype(radius) + pressure = similar(radius) + n_particles = length(radius) + + mass = ones(ELTYPE, n_particles) + + smoothing_length = 2 * maximum(radius) + + new{ELTYPE, Nothing, Nothing, Nothing, Nothing, Nothing}(pressure, mass, radius, nothing, nothing, nothing, smoothing_length, 1.0, nothing, nothing) end end From 4c67151f3262edee33c3a0d01f1832ab53958d9f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 26 Oct 2023 09:00:04 +0200 Subject: [PATCH 4/4] get basic simulation working --- examples/dem/rectangular_tank_2d.jl | 30 +++--- examples/fluid/rectangular_tank_2d.jl | 3 +- src/TrixiParticles.jl | 2 +- src/general/semidiscretization.jl | 14 ++- .../dummy_particles/dummy_particles.jl | 14 ++- src/schemes/boundary/rhs.jl | 2 +- src/schemes/boundary/system.jl | 67 +++++++++++-- .../solid/discrete_element_method/rhs.jl | 96 ++++++++++++++++++- .../solid/discrete_element_method/system.jl | 12 ++- src/visualization/write2vtk.jl | 2 +- 10 files changed, 201 insertions(+), 41 deletions(-) diff --git a/examples/dem/rectangular_tank_2d.jl b/examples/dem/rectangular_tank_2d.jl index b9df41024..5593db645 100644 --- a/examples/dem/rectangular_tank_2d.jl +++ b/examples/dem/rectangular_tank_2d.jl @@ -6,21 +6,21 @@ gravity = -9.81 # ========================================================================================== # ==== Fluid -particle_spacing = 0.25 +particle_spacing = 0.1 # Ratio of fluid particle spacing to boundary particle spacing beta = 1 -boundary_layers = 1 +boundary_layers = 2 -water_width = 2.0 -water_height = 2.0 -water_density = 1000.0 +rock_width = 2.0 +rock_height = 2.0 +rock_density = 3000.0 tank_width = 2.0 tank_height = 4.0 -tank = RectangularTank(particle_spacing, (water_width, water_height), - (tank_width, tank_height), water_density, +tank = RectangularTank(particle_spacing, (rock_width, rock_height), + (tank_width, tank_height), rock_density, n_layers=boundary_layers, spacing_ratio=beta) # ========================================================================================== @@ -31,16 +31,18 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.radius) # ========================================================================================== # ==== Systems -solid_system = DEMSystem(tank.fluid, 0.1, acceleration=(0.0, gravity)) - -# boundary_system = BoundaryDEMSystem(tank.boundary, boundary_model) +# let them fall +tank.fluid.coordinates[2,:] .+=0.5 +solid_system = DEMSystem(tank.fluid, 2*10^5, acceleration=(0.0, gravity)) +boundary_system = BoundaryDEMSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(solid_system, neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(solid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch) -tspan = (0.0, 2.0) +tspan = (0.0, 5.0) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) @@ -58,6 +60,6 @@ callbacks = CallbackSet(info_callback, saving_callback) # and the time integration method interprets this as an instability. sol = solve(ode, RDPK3SpFSAL49(), abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) - dtmax=1e-2, # Limit stepsize to prevent crashing + 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 save_everystep=false, callback=callbacks); diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl index 55f318ba9..87ae3b9da 100644 --- a/examples/fluid/rectangular_tank_2d.jl +++ b/examples/fluid/rectangular_tank_2d.jl @@ -40,7 +40,8 @@ tank = RectangularTank(particle_spacing, (water_width, water_height), boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, AdamiPressureExtrapolation(), - smoothing_kernel, particle_spacing, smoothing_factor) + smoothing_kernel, particle_spacing, + smoothing_factor) # K = 9.81 * water_height # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta, diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 8737c1e10..31faf8096 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -34,7 +34,7 @@ include("visualization/write2vtk.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, DEMSystem + BoundarySPHSystem, DEMSystem, BoundaryDEMSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 48006a3f4..7ac2944d6 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -106,6 +106,15 @@ end return compact_support(smoothing_kernel, smoothing_length) end +@inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem) + return 0 +end + +@inline function compact_support(system::BoundaryDEMSystem, neighbor::DEMSystem) + # use the compact support of the DEMSystem + return compact_support(neighbor, system) +end + @inline function compact_support(system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem) (; smoothing_kernel, smoothing_length) = system @@ -379,7 +388,10 @@ end return dv end -@inline add_acceleration!(dv, particle, system::BoundarySPHSystem) = dv +@inline function add_acceleration!(dv, particle, + system::Union{BoundarySPHSystem, BoundaryDEMSystem}) + dv +end @inline function add_damping_force!(dv, damping_coefficient, v, particle, system::FluidSystem) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 9d0aba397..164ba569a 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -78,7 +78,8 @@ struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C} cache :: C function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, - density_calculator, smoothing_kernel, particle_spacing, + density_calculator, smoothing_kernel, + particle_spacing, smoothing_factor; viscosity=NoViscosity(), state_equation=nothing) ELTYPE = eltype(initial_density) @@ -94,8 +95,10 @@ struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C} new{ELTYPE, typeof(state_equation), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(cache)}(pressure, hydrodynamic_mass, radius, state_equation, density_calculator, - smoothing_kernel, smoothing_length, smoothing_factor, viscosity, cache) + typeof(cache)}(pressure, hydrodynamic_mass, radius, state_equation, + density_calculator, + smoothing_kernel, smoothing_length, smoothing_factor, viscosity, + cache) end function BoundaryModelDummyParticles(radius) @@ -107,7 +110,10 @@ struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C} smoothing_length = 2 * maximum(radius) - new{ELTYPE, Nothing, Nothing, Nothing, Nothing, Nothing}(pressure, mass, radius, nothing, nothing, nothing, smoothing_length, 1.0, nothing, nothing) + new{ELTYPE, Nothing, Nothing, Nothing, Nothing, Nothing}(pressure, mass, radius, + nothing, nothing, nothing, + smoothing_length, 1.0, + nothing, nothing) end end diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index ddadfca4e..4bfd94f53 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -1,7 +1,7 @@ # Interaction of boundary with other systems function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::BoundarySPHSystem, + particle_system::Union{BoundarySPHSystem, BoundaryDEMSystem}, neighbor_system) # TODO Solids and moving boundaries should be considered in the continuity equation return dv diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 7c512d854..269b39972 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -27,6 +27,46 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: System{NDIMS} end end +""" + BoundaryDEMSystem(inititial_condition) + +System for boundaries modeled by boundary particles. +The interaction between fluid and boundary particles is specified by the boundary model. + +""" +struct BoundaryDEMSystem{BM, NDIMS, ELTYPE <: Real} <: System{NDIMS} + coordinates :: Array{ELTYPE, 2} + boundary_model :: BM + + function BoundaryDEMSystem(inititial_condition, model) + coordinates = inititial_condition.coordinates + NDIMS = size(coordinates, 1) + + return new{typeof(model), NDIMS, eltype(coordinates)}(coordinates, model) + end +end + +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 + +function Base.show(io::IO, ::MIME"text/plain", system::BoundaryDEMSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "BoundaryDEMSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "boundary model", system.boundary_model) + summary_footer(io) + end +end + """ BoundaryMovement(movement_function, is_moving) @@ -90,13 +130,17 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) end end -timer_name(::BoundarySPHSystem) = "boundary" +timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary" -@inline Base.eltype(system::BoundarySPHSystem) = eltype(system.coordinates) +@inline function Base.eltype(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) + eltype(system.coordinates) +end # This does not account for moving boundaries, but it's only used to initialize the # neighborhood search, anyway. -@inline initial_coordinates(system::BoundarySPHSystem) = system.coordinates +@inline function initial_coordinates(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) + system.coordinates +end function (movement::BoundaryMovement)(system, t) (; coordinates, cache) = system @@ -128,13 +172,13 @@ function (movement::Nothing)(system, t) return system end -@inline function nparticles(system::BoundarySPHSystem) - length(system.boundary_model.hydrodynamic_mass) +@inline function nparticles(system::Union{BoundaryDEMSystem, BoundarySPHSystem}) + size(system.coordinates, 2) end # No particle positions are advanced for boundary systems, # except when using BoundaryModelDummyParticles with ContinuityDensity. -@inline function n_moving_particles(system::BoundarySPHSystem) +@inline function n_moving_particles(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) return 0 end @@ -152,13 +196,14 @@ end nparticles(system) end -@inline u_nvariables(system::BoundarySPHSystem) = 0 +@inline u_nvariables(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) = 0 # For BoundaryModelDummyParticles with ContinuityDensity, this needs to be 1. # For all other models and density calculators, it's irrelevant. @inline v_nvariables(system::BoundarySPHSystem) = 1 +@inline v_nvariables(system::BoundaryDEMSystem) = 0 -@inline function current_coordinates(u, system::BoundarySPHSystem) +@inline function current_coordinates(u, system::Union{BoundarySPHSystem, BoundaryDEMSystem}) return system.coordinates end @@ -237,11 +282,13 @@ function update_final!(system::BoundarySPHSystem, system_index, v, u, v_ode, u_o return system end -function write_u0!(u0, system::BoundarySPHSystem) +function write_u0!(u0, system::Union{BoundarySPHSystem, BoundaryDEMSystem}) return u0 end -function write_v0!(v0, system::BoundarySPHSystem{<:BoundaryModelMonaghanKajtar}) +function write_v0!(v0, + system::Union{BoundarySPHSystem{<:BoundaryModelMonaghanKajtar}, + BoundaryDEMSystem}) return v0 end diff --git a/src/schemes/solid/discrete_element_method/rhs.jl b/src/schemes/solid/discrete_element_method/rhs.jl index b6464c625..d2b2b3daf 100644 --- a/src/schemes/solid/discrete_element_method/rhs.jl +++ b/src/schemes/solid/discrete_element_method/rhs.jl @@ -7,8 +7,21 @@ # The magnitude of the force is determined by the stiffness constant `kn` and the overlap distan function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::DEMSystem, - neighbor_system::DEMSystem) - (; mass, radius, kn) = particle_system + neighbor_system) + (; mass, radius) = particle_system + nghb_radius = neighbor_system.radius + nghb_mass = neighbor_system.mass + gamma_coefficient = 0.0001 + + E_particle = 10*10^9 #particle_system.elastic_modulus # Elastic modulus for particles + nu_particle = 0.3 #particle_system.poissons_ratio # Poisson's ratio for particles + + # Extracting material properties for neighbor system + E_neighbor = 10*10^9 #neighbor_system.elastic_modulus # Elastic modulus for neighbor + nu_neighbor = 0.3 #neighbor_system.poissons_ratio # Poisson's ratio for neighbor + + # Compute effective modulus for both systems + E_star = 1 / ((1 - nu_particle^2) / E_particle + (1 - nu_neighbor^2) / E_neighbor) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -19,16 +32,35 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, distance < sqrt(eps()) && return # Calculate the overlap (penetration depth) between the two particles - overlap = radius[particle] + radius[neighbor] - distance + overlap = radius[particle] + nghb_radius[neighbor] - distance # If there's no overlap, no force needs to be applied overlap <= 0 && return + # Compute effective radius for the interaction + r_star = (radius[particle] * nghb_radius[neighbor]) / (radius[particle] + nghb_radius[neighbor]) + + # Compute stiffness constant kn for the interaction + kn = (4/3) * E_star * sqrt(r_star * overlap) + + # Calculate effective mass for the interaction + m_star = (mass[particle] * nghb_mass[neighbor]) / (mass[particle] + nghb_mass[neighbor]) + + # Calculate critical damping coefficient + gamma_c = 2 * sqrt(m_star * kn) + # Normal direction from particle to neighbor normal = pos_diff / distance - # Compute the force magnitude using Hertzian contact mechanics - force_magnitude = kn * overlap + part_v = extract_svector(v_particle_system, particle_system, particle) + nghbr_v = extract_svector(v_neighbor_system, neighbor_system, neighbor) + + # Relative velocity + rel_vel = part_v - nghbr_v + rel_vel_normal = dot(rel_vel, normal) + + # Compute the force magnitude using Hertzian contact mechanics with damping + force_magnitude = kn * overlap + gamma_coefficient * gamma_c * rel_vel_normal force = force_magnitude * normal @@ -40,3 +72,57 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, return dv end + +function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, + u_neighbor_system, neighborhood_search, particle_system::DEMSystem, + neighbor_system::BoundaryDEMSystem) + (; mass, radius, kn) = particle_system + + nghb_radius = neighbor_system.boundary_model.radius + + max_kn = kn + wall_kn_factor = 5 * max_kn + wall_kn = wall_kn_factor * max_kn + + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + for_particle_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + # Calculate the overlap (penetration depth) between the two particles + overlap = radius[particle] + nghb_radius[neighbor] - distance + + # If there's no overlap, no force needs to be applied + overlap <= 0 && return + + # Normal direction from particle to neighbor + normal = pos_diff / distance + + # Position correction to prevent penetration + position_correction_factor = 0.5 # you can tweak this value + @inbounds for i in 1:ndims(particle_system) + u_particle_system[i, particle] -= position_correction_factor * overlap * normal[i] + end + + # Compute the force magnitude using Hertzian contact mechanics + force_magnitude = wall_kn * position_correction_factor * overlap + + force = force_magnitude * normal + + # Reset the velocity component in the normal direction if it's pushing the particle into the wall + normal_velocity = dot(dv[:, particle], normal) + if normal_velocity < 0 + @inbounds for i in 1:ndims(particle_system) + dv[i, particle] -= normal_velocity * normal[i] + # we still need to push against the particles coming to close + dv[i, particle] += force[i] / mass[particle] + end + end + end + + return dv +end diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index cc666cf3d..dc252ba2b 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -5,8 +5,9 @@ struct DEMSystem{NDIMS, ELTYPE <: Real} <: System{NDIMS} kn :: ELTYPE # Normal stiffness acceleration :: SVector{NDIMS, ELTYPE} - function DEMSystem(initial_condition, kn; acceleration=ntuple(_ -> 0.0, - ndims(initial_condition))) + function DEMSystem(initial_condition, kn; + acceleration=ntuple(_ -> 0.0, + ndims(initial_condition))) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -19,7 +20,6 @@ struct DEMSystem{NDIMS, ELTYPE <: Real} <: System{NDIMS} throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - return new{NDIMS, ELTYPE}(initial_condition, mass, radius, kn, acceleration_) end end @@ -66,3 +66,9 @@ function compact_support(system::DEMSystem, neighbor::DEMSystem) # todo: needs to be changed for more complex simulations return 3 * max(maximum(system.radius), maximum(neighbor.radius)) end + +function compact_support(system::DEMSystem, neighbor) + # we for now assume that the compact support is 3 * radius + # todo: needs to be changed for more complex simulations + return 3 * maximum(system.radius) +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index e9dcb9302..325761780 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -140,7 +140,7 @@ end vtkname(system) = "particle" vtkname(system::FluidSystem) = "fluid" vtkname(system::Union{TotalLagrangianSPHSystem, DEMSystem}) = "solid" -vtkname(system::BoundarySPHSystem) = "boundary" +vtkname(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary" function write2vtk!(vtk, v, u, t, system; write_meta_data=true) vtk["velocity"] = view(v, 1:ndims(system), :)