diff --git a/NEWS.md b/NEWS.md index 3541c9d68..ffe272e8a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -17,17 +17,16 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between ## Version 0.1.2 +### Added +A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH + +## Version 0.1.1 + ### Highlights #### Discrete Element Method A basic implementation of the discrete element method was added. - -## Version 0.1.1 - -### Added -A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids", 2013 was added to WCSPH - # Pre Initial Release (v0.1.0) This section summarizes the initial features that TrixiParticles.jl was released with. diff --git a/Project.toml b/Project.toml index 20b40306c..ea1bc9aa5 100644 --- a/Project.toml +++ b/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.1.2-dev" +version = "0.1.3-dev" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/src/index.md b/docs/src/index.md index a3b41a417..42c2c42c6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,6 +5,7 @@ TrixiParticles.jl is a numerical simulation framework designed for particle-base ## Features - Incompressible Navier-Stokes - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC) + - Models: Surface Tension - Solid-body mechanics - Methods: Total Lagrangian SPH (TLSPH) - Fluid-Structure Interaction diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 69b514950..353b6b5cb 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -147,3 +147,65 @@ Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "density_diffus Modules = [TrixiParticles] Pages = [joinpath("general", "corrections.jl")] ``` + +## [Surface Tension](@id surface_tension) + +### Akinci-based intra-particle force surface tension and wall adhesion model +The work by Akinci proposes three forces: +- a cohesion force +- a surface area minimization force +- a wall adhesion force + +The classical model is composed of the curvature minimization and cohesion force. + +#### Cohesion force +The model calculates the cohesion force based on the distance between particles and the support radius ``h_c``. +This force is determined using two distinct regimes within the support radius: +- For particles closer than half the support radius, + a repulsive force is calculated to prevent particle clustering too tightly, + enhancing the simulation's stability and realism. +- Beyond half the support radius and within the full support radius, + an attractive force is computed, simulating the effects of surface tension that draw particles together. +The cohesion force, ``F_{\text{cohesion}}``, for a pair of particles is given by: +```math +F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, +``` +where: +- ``\sigma`` represents the surface tension coefficient, adjusting the overall strength of the cohesion effect. +- ``C`` is a scalar function of the distance between particles. + +The cohesion kernel ``C`` is defined as +```math +C(r)=\frac{32}{\pi h_c^9} +\begin{cases} +(h_c-r)^3 r^3, & \text{if } 2r > h_c \\ +2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c \\ +0, & \text{otherwise} +\end{cases} +``` + +#### Surface area minimization force +To model the minimization of the surface area and curvature of the fluid, a curvature force is used, which is calculated as +```math +F_{\text{curvature}} = -\sigma (n_a - n_b) +``` + +#### Wall adhesion force +The wall adhesion model proposed by Akinci et al. is based on a kernel function which is 0 from 0.0 to 0.5 support radiia with a maximum at 0.75. +With the force calculated with an adhesion coefficient ``\beta`` as +```math +F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, +``` +with ``A`` being the adhesion kernel defined as +```math +A(r)= \frac{0.007}{h_c^{3.25}} +\begin{cases} +\sqrt[4]{- \frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c \\ +0, & \text{otherwise.} +\end{cases} +``` + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "surface_tension.jl")] +``` diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 40176927c..e58a7cc06 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -57,8 +57,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, density_diffusion=density_diffusion, - acceleration=(0.0, -gravity), - correction=nothing) + acceleration=(0.0, -gravity), correction=nothing, + surface_tension=nothing) # ========================================================================================== # ==== Boundary @@ -69,7 +69,7 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar smoothing_kernel, smoothing_length, correction=nothing) -boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0) # ========================================================================================== # ==== Simulation diff --git a/examples/fluid/dam_break_2d_surface_tension.jl b/examples/fluid/dam_break_2d_surface_tension.jl new file mode 100644 index 000000000..d83eb9865 --- /dev/null +++ b/examples/fluid/dam_break_2d_surface_tension.jl @@ -0,0 +1,31 @@ +# This example shows how surface tension can be applied to existing cases +# and which parameters have to be changed! +using TrixiParticles + +fluid_density = 1000.0 + +H = 0.6 +fluid_particle_spacing = H / 60 + +# Set the surface tension to a value that is accurate in your case. +# Note: This usually requires calibration to be physically accurate! +surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.025) + +# `density_diffusion` is deactivated since the interaction with the surface tension model can +# cause stability problems. +# `adhesion_coefficient` needs to be set to a value so that the fluid doesn't separate +# from the boundary. +# Note: The surface tension model leads to an increase in compressibility of the fluid, +# which needs to be rectified by an increase of the `sound_speed`. +# Note: The Wendland Kernels don't work very well here since the `SurfaceTensionAkinci` +# model is optimized for a compact support of `2 * particle_spacing`, which would result +# in a too small `smoothing_length` for the Wendland Kernel functions. +# Note: Adhesion will result in additional friction at the boundary. +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + surface_tension=surface_tension, + fluid_particle_spacing=fluid_particle_spacing, + smoothing_kernel=SchoenbergCubicSplineKernel{2}(), + smoothing_length=1.0 * fluid_particle_spacing, + correction=AkinciFreeSurfaceCorrection(fluid_density), + density_diffusion=nothing, adhesion_coefficient=0.05, + sound_speed=100.0, tspan=(0.0, 2.0)) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl new file mode 100644 index 000000000..37132d70b --- /dev/null +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -0,0 +1,95 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.3) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.2) +sphere2_center = (1.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 1.0 * fluid_particle_spacing +fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.005 +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, + viscosity=viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), + correction=AkinciFreeSurfaceCorrection(fluid_density)) + +sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=0.001) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="", + write_meta_data=true) + +callbacks = CallbackSet(info_callback, saving_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl new file mode 100644 index 000000000..a4ebb799e --- /dev/null +++ b/examples/fluid/falling_water_spheres_3d.jl @@ -0,0 +1,37 @@ +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.008 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +nu = 0.01 +fluid_density = 1000.0 +sound_speed = 50 + +sphere1_radius = 0.05 + +sphere1_center = (0.5, 0.5, 0.2) +sphere2_center = (1.5, 0.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0)) + +# `compact_support` needs to be `2.0 * particle_spacing` to be correct +fluid_smoothing_length = 1.0 * fluid_particle_spacing + +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), + fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.2), + initial_fluid_size=(0.0, 0.0, 0.0), + tank_size=(2.0, 1.0, 0.1), sound_speed=sound_speed, + faces=(true, true, true, true, true, false), + acceleration=(0.0, 0.0, -gravity), sphere1=sphere1, sphere2=sphere2, + fluid_smoothing_length=fluid_smoothing_length, + fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(), + nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed), + surface_tension_coefficient=1.5, adhesion_coefficient=0.5) diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl new file mode 100644 index 000000000..9006ebf0e --- /dev/null +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -0,0 +1,59 @@ +# In this example we can observe that the `SurfaceTensionAkinci` surface tension model correctly leads to a +# surface minimization of the water square and approaches a sphere. +using TrixiParticles +using OrdinaryDiffEq + +fluid_density = 1000.0 + +particle_spacing = 0.1 + +# Note: Only square shapes will result in a sphere. +# Furthermore, changes of the coefficients might be necessary for higher resolutions or larger squares. +fluid_size = (0.5, 0.5) + +sound_speed = 20.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=true) + +# For all surface tension simulations, we need a compact support of `2 * particle_spacing` +# smoothing_length = 2.0 * particle_spacing +# smoothing_kernel = WendlandC2Kernel{2}() +# nu = 0.01 + +smoothing_length = 1.0 * particle_spacing +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +nu = 0.025 + +fluid = RectangularShape(particle_spacing, round.(Int, fluid_size ./ particle_spacing), + zeros(length(fluid_size)), density=fluid_density) + +alpha = 8 * nu / (smoothing_length * sound_speed) +source_terms = SourceTermDamping(; damping_coefficient=0.5) +fluid_system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), + state_equation, smoothing_kernel, + smoothing_length, + viscosity=ArtificialViscosityMonaghan(alpha=alpha, + beta=0.0), + surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02), + correction=AkinciFreeSurfaceCorrection(fluid_density), + source_terms=source_terms) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system) + +tspan = (0.0, 3.0) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +# For overwriting via `trixi_include` +saving_callback = SolutionSavingCallback(dt=0.02) + +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); diff --git a/examples/fluid/sphere_surface_tension_3d.jl b/examples/fluid/sphere_surface_tension_3d.jl new file mode 100644 index 000000000..58a506294 --- /dev/null +++ b/examples/fluid/sphere_surface_tension_3d.jl @@ -0,0 +1,25 @@ +# In this example we can observe that the `SurfaceTensionAkinci` surface tension model correctly leads to a +# surface minimization of the water cube and approaches a sphere. +using TrixiParticles +using OrdinaryDiffEq + +fluid_density = 1000.0 + +particle_spacing = 0.1 +fluid_size = (0.9, 0.9, 0.9) + +sound_speed = 20.0 + +# For all surface tension simulations, we need a compact support of `2 * particle_spacing` +smoothing_length = 1.0 * particle_spacing + +nu = 0.01 + +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "sphere_surface_tension_2d.jl"), + dt=0.1, cfl=1.2, surface_tension_coefficient=0.1, + tspan=(0.0, 10.0), nu=nu, + alpha=10 * nu / (smoothing_length * sound_speed), + smoothing_kernel=SchoenbergCubicSplineKernel{3}(), + particle_spacing=particle_spacing, sound_speed=sound_speed, + fluid_density=fluid_density, fluid_size=fluid_size) diff --git a/examples/fluid/sphere_surface_tension_wall_2d.jl b/examples/fluid/sphere_surface_tension_wall_2d.jl new file mode 100644 index 000000000..3e1cbdb80 --- /dev/null +++ b/examples/fluid/sphere_surface_tension_wall_2d.jl @@ -0,0 +1,44 @@ +# In this example we try to approach the static shape of a water droplet on a horizontal plane. +# The shape of a static droplet can be calculated from the Young-Laplace equation. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.0025 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 0.5) + +# Boundary geometry and initial fluid particle positions +tank_size = (0.5, 0.1) + +fluid_density = 1000.0 +sound_speed = 120.0 + +sphere_radius = 0.05 + +sphere1_center = (0.25, sphere_radius) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere()) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 1.0 * fluid_particle_spacing + +# For perfect wetting +# nu = 0.0005 +# For no wetting +nu = 0.001 + +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +# `adhesion_coefficient = 1.0` and `surface_tension_coefficient = 0.01` for perfect wetting +# `adhesion_coefficient = 0.001` and `surface_tension_coefficient = 2.0` for no wetting +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"), + sphere=nothing, sphere1=sphere1, adhesion_coefficient=0.001, + wall_viscosity=4.0 * nu, surface_tension_coefficient=2.0, alpha=alpha, + sound_speed=sound_speed, fluid_density=fluid_density, nu=nu, + fluid_particle_spacing=fluid_particle_spacing, tspan=tspan, + tank_size=tank_size) diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index 17868d201..04c1f1af6 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -9,8 +9,9 @@ using Polyester include("n_body_system.jl") -# Redefine interact in a more optimized way. -function TrixiParticles.interact!(du, u_particle_system, u_neighbor_system, +# Redefine interact in a more optimized way +function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::NBodySystem, neighbor_system::NBodySystem) @@ -34,16 +35,15 @@ function TrixiParticles.interact!(du, u_particle_system, u_neighbor_system, tmp1 = mass[neighbor] * tmp tmp2 = mass[particle] * tmp - for i in 1:ndims(particle_system) - j = i + ndims(particle_system) - # This is slightly faster than du[...] += tmp1 * pos_diff[i] - du[j, particle] = muladd(tmp1, pos_diff[i], du[j, particle]) - du[j, neighbor] = muladd(tmp2, -pos_diff[i], du[j, neighbor]) + @inbounds for i in 1:ndims(particle_system) + # This is slightly faster than dv[...] += tmp1 * pos_diff[i] + dv[i, particle] = muladd(tmp1, pos_diff[i], dv[i, particle]) + dv[i, neighbor] = muladd(tmp2, -pos_diff[i], dv[i, neighbor]) end end end - return du + return dv end # ========================================================================================== @@ -98,6 +98,8 @@ function symplectic_euler!(velocity, coordinates, semi) u[i] += 0.01 * du[i] end end + + return v, u end # One RHS evaluation is so fast that timers make it multiple times slower. diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 96b6cdebf..caf40a62e 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -1,7 +1,7 @@ using TrixiParticles using LinearAlgebra -struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS} +struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothing} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] G :: ELTYPE diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c95451a25..f170e08d8 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -70,5 +70,6 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk +export SurfaceTensionAkinci, CohesionForceAkinci end # module diff --git a/src/general/general.jl b/src/general/general.jl index 4c9531683..dfe52e06a 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -1,12 +1,17 @@ -abstract type System{NDIMS} end +# Abstract supertype for all system types. We additionally store the type of the system's +# initial condition, which is `Nothing` when using KernelAbstractions.jl. +abstract type System{NDIMS, IC} end -abstract type FluidSystem{NDIMS} <: System{NDIMS} end +# When using KernelAbstractions.jl, the initial condition has been replaced by `nothing` +GPUSystem = System{NDIMS, Nothing} where {NDIMS} + +abstract type FluidSystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::FluidSystem) = "fluid" -abstract type SolidSystem{NDIMS} <: System{NDIMS} end +abstract type SolidSystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::SolidSystem) = "solid" -abstract type BoundarySystem{NDIMS} <: System{NDIMS} end +abstract type BoundarySystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::BoundarySystem) = "boundary" @inline function set_zero!(du) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a3ce1e056..6a1d703b2 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -60,6 +60,8 @@ struct Semidiscretization{S, RU, RV, NS} end end +GPUSemidiscretization = Semidiscretization{<:NTuple{<:Any, GPUSystem}} + function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch, periodic_box_min_corner=nothing, periodic_box_max_corner=nothing, threaded_nhs_update=true) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 843dfa4df..cade5fbe9 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -1,5 +1,5 @@ """ - BoundarySPHSystem(initial_condition, boundary_model; movement=nothing) + BoundarySPHSystem(initial_condition, boundary_model; movement=nothing, adhesion_coefficient=0.0) System for boundaries modeled by boundary particles. The interaction between fluid and boundary particles is specified by the boundary model. @@ -10,30 +10,37 @@ The interaction between fluid and boundary particles is specified by the boundar # Keyword Arguments - `movement`: For moving boundaries, a [`BoundaryMovement`](@ref) can be passed. +- `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface. + Note: currently it is assumed that all fluids have the same adhesion coefficient. """ -struct BoundarySPHSystem{BM, NDIMS, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS} - initial_condition :: IC - coordinates :: CO # Array{ELTYPE, 2} - boundary_model :: BM - movement :: M - ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) - cache :: CA +struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS, IC} + initial_condition :: IC + coordinates :: CO # Array{ELTYPE, 2} + boundary_model :: BM + movement :: M + ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) + adhesion_coefficient :: ELTYPE + cache :: CA # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. function BoundarySPHSystem(initial_condition, coordinates, boundary_model, movement, - ismoving, cache) - new{typeof(boundary_model), size(coordinates, 1), + ismoving, adhesion_coefficient, cache) + ELTYPE = eltype(coordinates) + + new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), typeof(coordinates), typeof(movement), typeof(ismoving), typeof(cache)}(initial_condition, coordinates, boundary_model, - movement, - ismoving, cache) + movement, ismoving, + adhesion_coefficient, cache) end end -function BoundarySPHSystem(initial_condition, model; movement=nothing) +function BoundarySPHSystem(initial_condition, model; movement=nothing, + adhesion_coefficient=0.0) coordinates = copy(initial_condition.coordinates) + ismoving = Ref(!isnothing(movement)) cache = create_cache_boundary(movement, initial_condition) @@ -45,8 +52,9 @@ function BoundarySPHSystem(initial_condition, model; movement=nothing) movement.moving_particles .= collect(1:nparticles(initial_condition)) end - return BoundarySPHSystem(initial_condition, coordinates, model, - movement, ismoving, cache) + # Because of dispatches boundary model needs to be first! + return BoundarySPHSystem(initial_condition, coordinates, model, movement, + ismoving, adhesion_coefficient, cache) end """ @@ -59,10 +67,12 @@ The interaction between fluid and boundary particles is specified by the boundar This is an experimental feature and may change in a future releases. """ -struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D} <: BoundarySystem{NDIMS} - coordinates :: ARRAY2D # [dimension, particle] - radius :: ARRAY1D # [particle] - normal_stiffness :: ELTYPE +struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC, + ARRAY1D, ARRAY2D} <: BoundarySystem{NDIMS, IC} + initial_condition :: IC + coordinates :: ARRAY2D # [dimension, particle] + radius :: ARRAY1D # [particle] + normal_stiffness :: ELTYPE function BoundaryDEMSystem(initial_condition, normal_stiffness) coordinates = initial_condition.coordinates @@ -70,9 +80,9 @@ struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D} <: BoundarySys ones(length(initial_condition.mass)) NDIMS = size(coordinates, 1) - return new{NDIMS, eltype(coordinates), typeof(radius), typeof(coordinates)}(coordinates, - radius, - normal_stiffness) + return new{NDIMS, eltype(coordinates), typeof(initial_condition), + typeof(radius), typeof(coordinates)}(initial_condition, coordinates, + radius, normal_stiffness) end end @@ -156,6 +166,7 @@ function Base.show(io::IO, system::BoundarySPHSystem) print(io, "BoundarySPHSystem{", ndims(system), "}(") print(io, system.boundary_model) print(io, ", ", system.movement) + print(io, ", ", system.adhesion_coefficient) print(io, ") with ", nparticles(system), " particles") end @@ -171,6 +182,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem) summary_line(io, "movement function", isnothing(system.movement) ? "nothing" : string(system.movement.movement_function)) + summary_line(io, "adhesion coefficient", system.adhesion_coefficient) summary_footer(io) end end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 15c765040..f59bca7d3 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -40,7 +40,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, - PF, ST, C} <: FluidSystem{NDIMS} + PF, ST, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 04e29888a..056662bbc 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -70,5 +70,6 @@ end include("pressure_acceleration.jl") include("viscosity.jl") +include("surface_tension.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl new file mode 100644 index 000000000..69e468cd7 --- /dev/null +++ b/src/schemes/fluid/surface_tension.jl @@ -0,0 +1,206 @@ +abstract type AkinciTypeSurfaceTension end + +@doc raw""" + CohesionForceAkinci(surface_tension_coefficient=1.0) + +This model only implements the cohesion force of the Akinci surface tension model. + +# Keywords +- `surface_tension_coefficient=1.0`: Modifies the intensity of the surface tension-induced force, + enabling the tuning of the fluid's surface tension properties within the simulation. + +# References +- Nadir Akinci, Gizem Akinci, Matthias Teschner. + "Versatile Surface Tension and Adhesion for SPH Fluids". + In: ACM Transactions on Graphics 32.6 (2013). + [doi: 10.1145/2508363.2508395](https://doi.org/10.1145/2508363.2508395) +""" +struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension + surface_tension_coefficient::ELTYPE + + function CohesionForceAkinci(; surface_tension_coefficient=1.0) + new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + end +end + +@doc raw""" + SurfaceTensionAkinci(surface_tension_coefficient=1.0) + +Implements a model for surface tension and adhesion effects drawing upon the +principles outlined by Akinci et al. This model is instrumental in capturing the nuanced +behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or +separation, by utilizing intra-particle forces. + +# Keywords +- `surface_tension_coefficient=1.0`: A parameter to adjust the magnitude of + surface tension forces, facilitating the fine-tuning of how surface tension phenomena + are represented in the simulation. + +# References +- Nadir Akinci, Gizem Akinci, Matthias Teschner. + "Versatile Surface Tension and Adhesion for SPH Fluids". + In: ACM Transactions on Graphics 32.6 (2013). + [doi: 10.1145/2508363.2508395](https://doi.org/10.1145/2508363.2508395) +""" +struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension + surface_tension_coefficient::ELTYPE + + function SurfaceTensionAkinci(; surface_tension_coefficient=1.0) + new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + end +end + +# Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`. +# Currently, specializations reducing this to simple multiplications exist only up +# to a power of three, see +# https://github.com/JuliaLang/julia/blob/34934736fa4dcb30697ac1b23d11d5ad394d6a4d/base/intfuncs.jl#L327-L339 +# By using the `@fastpow` macro, we are consciously trading off some precision in the result +# for enhanced computational speed. This is especially useful in scenarios where performance +# is a higher priority than exact precision. +@fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, distance) + (; surface_tension_coefficient) = surface_tension + + # Eq. 2 + # We only reach this function when `sqrt(eps()) < distance <= support_radius` + if distance > 0.5 * support_radius + # Attractive force + C = (support_radius - distance)^3 * distance^3 + else + # `distance < 0.5 * support_radius` + # Repulsive force + C = 2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64.0 + end + C *= 32.0 / (pi * support_radius^9) + + # Eq. 1 in acceleration form + cohesion_force = -surface_tension_coefficient * m_b * C * pos_diff / distance + + return cohesion_force +end + +@inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, + distance, adhesion_coefficient) + + # The neighborhood search has an `<=` check, but for `distance == support_radius` + # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`. + # TODO Change this in the neighborhood search? + # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19 + distance >= support_radius && return zero(pos_diff) + + distance <= 0.5 * support_radius && return zero(pos_diff) + + # Eq. 7 + A = 0.007 / support_radius^3.25 * + (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25 + + # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as + # `m_b = rho_0 * volume` (Akinci boundary condition treatment) + adhesion_force = -adhesion_coefficient * m_b * A * pos_diff / distance + + return adhesion_force +end + +# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" +# Note: most of the time this only leads to an approximation of the surface normal + +function calc_normal_akinci!(system, neighbor_system::FluidSystem, + surface_tension::SurfaceTensionAkinci, u_system, + v_neighbor_system, u_neighbor_system, + neighborhood_search) + (; smoothing_length, cache) = system + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + for_particle_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + density_neighbor = particle_density(v_neighbor_system, + neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, + particle) + for i in 1:ndims(system) + cache.surface_normal[i, particle] += m_b / density_neighbor * + grad_kernel[i] * smoothing_length + end + end + + return system +end + +function calc_normal_akinci!(system, neighbor_system, surface_tension, u_system, + v_neighbor_system, u_neighbor_system, + neighborhood_search) + # Normal not needed + return system +end + +@inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, + surface_tension_b::CohesionForceAkinci, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance) + (; smoothing_length) = particle_system + # No cohesion with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + support_radius = compact_support(smoothing_kernel, smoothing_length) + + return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) +end + +@inline function surface_tension_force(surface_tension_a::SurfaceTensionAkinci, + surface_tension_b::SurfaceTensionAkinci, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance) + (; smoothing_length, smoothing_kernel) = particle_system + (; surface_tension_coefficient) = surface_tension_a + + # No surface tension with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + n_a = surface_normal(surface_tension_a, particle_system, particle) + n_b = surface_normal(surface_tension_b, neighbor_system, neighbor) + support_radius = compact_support(smoothing_kernel, smoothing_length) + + return cohesion_force_akinci(surface_tension_a, support_radius, m_b, + pos_diff, distance) .- + (surface_tension_coefficient * (n_a - n_b)) +end + +# Skip +@inline function surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, particle, neighbor, + pos_diff, distance) + return zero(pos_diff) +end + +@inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, + particle_system::FluidSystem, + neighbor_system::BoundarySystem, particle, neighbor, + pos_diff, distance) + (; smoothing_length, smoothing_kernel) = particle_system + (; adhesion_coefficient, boundary_model) = neighbor_system + + # No adhesion with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + # No reason to calculate the adhesion force if adhesion coefficient is near zero + abs(adhesion_coefficient) < eps() && return zero(pos_diff) + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + support_radius = compact_support(smoothing_kernel, smoothing_length) + return adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, + adhesion_coefficient) +end + +@inline function adhesion_force(surface_tension, particle_system, neighbor_system, particle, + neighbor, pos_diff, distance) + return zero(pos_diff) +end diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 6a34692f0..554f41c14 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -190,6 +190,7 @@ end rho_a = particle_density(v_particle_system, particle_system, particle) rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + # TODO This is not correct for two different fluids. It should be `nu_a` and `nu_b`. eta_a = nu * rho_a eta_b = nu * rho_b diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index a1710052d..716245be8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -224,3 +224,35 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search return density_diffusion end + +@inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::FluidSystem, + neighbor_system::FluidSystem, + grad_kernel) + # Density diffusion terms are all zero for distance zero + distance < sqrt(eps()) && return + + (; delta) = density_diffusion + (; smoothing_length, state_equation) = particle_system + (; sound_speed) = state_equation + + volume_b = m_b / rho_b + + psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance, + particle_system, particle, neighbor) + density_diffusion_term = dot(psi, grad_kernel) * volume_b + + dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term +end + +# Density diffusion `nothing` or interaction other than fluid-fluid +@inline function density_diffusion!(dv, density_diffusion, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system, neighbor_system, grad_kernel) + return dv +end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 0b913bc23..5f1c41912 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,9 +6,12 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, particle_system::WeaklyCompressibleSPHSystem, neighbor_system) - (; density_calculator, state_equation, correction) = particle_system + (; density_calculator, state_equation, correction, surface_tension) = particle_system (; sound_speed) = state_equation + surface_tension_a = surface_tension_model(particle_system) + surface_tension_b = surface_tension_model(neighbor_system) + system_coords = current_coordinates(u_particle_system, particle_system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -25,9 +28,9 @@ function interact!(dv, v_particle_system, u_particle_system, rho_mean = 0.5 * (rho_a + rho_b) # Determine correction values - viscosity_correction, pressure_correction = free_surface_correction(correction, - particle_system, - rho_mean) + viscosity_correction, pressure_correction, surface_tension_correction = free_surface_correction(correction, + particle_system, + rho_mean) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) @@ -55,8 +58,17 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_mean) - for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension = surface_tension_correction * + surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance) + + dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, + particle, neighbor, pos_diff, distance) + + @inbounds for i in 1:ndims(particle_system) + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + + dv_adhesion[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end @@ -104,38 +116,6 @@ end particle_system, neighbor_system, grad_kernel) end -@inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::WeaklyCompressibleSPHSystem, - neighbor_system::WeaklyCompressibleSPHSystem, - grad_kernel) - # Density diffusion terms are all zero for distance zero - distance < sqrt(eps()) && return - - (; delta) = density_diffusion - (; smoothing_length, state_equation) = particle_system - (; sound_speed) = state_equation - - volume_b = m_b / rho_b - - psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance, - particle_system, particle, neighbor) - density_diffusion_term = dot(psi, grad_kernel) * volume_b - - dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term -end - -# Density diffusion `nothing` or interaction other than fluid-fluid -@inline function density_diffusion!(dv, density_diffusion, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) - return dv -end - @inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system, particle_system, neighbor_system, particle, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 7df899868..cb9a1da87 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -37,10 +37,12 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. [`BoundaryModelDummyParticles`](@ref) and [`AdamiPressureExtrapolation`](@ref). The keyword argument `acceleration` should be used instead for gravity-like source terms. +- `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) + """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, C} <: FluidSystem{NDIMS} + V, DD, COR, PF, ST, SRFT, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -54,6 +56,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, correction :: COR pressure_acceleration_formulation :: PF source_terms :: ST + surface_tension :: SRFT cache :: C end @@ -66,7 +69,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - correction=nothing, source_terms=nothing) + correction=nothing, source_terms=nothing, + surface_tension=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -98,14 +102,16 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(correction, initial_condition.density, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., + cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, - acceleration_, viscosity, - density_diffusion, correction, - pressure_acceleration, - source_terms, cache) + smoothing_kernel, smoothing_length, acceleration_, + viscosity, density_diffusion, correction, + pressure_acceleration, source_terms, surface_tension, + cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -132,6 +138,11 @@ function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_p return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end +function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) + surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + return (; surface_normal) +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time @@ -142,6 +153,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) + print(io, ", ", system.surface_tension) print(io, ", ", system.acceleration) print(io, ", ", system.source_terms) print(io, ") with ", nparticles(system), " particles") @@ -163,6 +175,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity) summary_line(io, "density diffusion", system.density_diffusion) + summary_line(io, "surface tension", system.surface_tension) summary_line(io, "acceleration", system.acceleration) summary_line(io, "source terms", system.source_terms |> typeof |> nameof) summary_footer(io) @@ -201,7 +214,7 @@ function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction) = system + (; density_calculator, correction, surface_tension) = system compute_correction_values!(system, correction, u, v_ode, u_ode, semi) @@ -211,6 +224,7 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) compute_pressure!(system, v) + compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) return system end @@ -325,3 +339,39 @@ end @inline function correction_matrix(system::WeaklyCompressibleSPHSystem, particle) extract_smatrix(system.cache.correction_matrix, system, particle) end + +function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_surface_normal!(system, surface_tension::SurfaceTensionAkinci, v, u, v_ode, + u_ode, semi, t) + (; cache) = system + + # Reset surface normal + set_zero!(cache.surface_normal) + + @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + nhs = get_neighborhood_search(system, semi) + + calc_normal_akinci!(system, neighbor_system, surface_tension, u, v_neighbor_system, + u_neighbor_system, nhs) + end + return system +end + +@inline function surface_normal(::SurfaceTensionAkinci, particle_system::FluidSystem, + particle) + (; cache) = particle_system + return extract_svector(cache.surface_normal, particle_system, particle) +end + +@inline function surface_tension_model(system::WeaklyCompressibleSPHSystem) + return system.surface_tension +end + +@inline function surface_tension_model(system) + return nothing +end diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index a7cd2868f..5aadfcfcc 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -26,8 +26,8 @@ specified material properties and contact mechanics. !!! warning "Experimental Implementation" This is an experimental feature and may change in a future releases. """ -struct DEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ST} <: SolidSystem{NDIMS} - initial_condition :: InitialCondition{ELTYPE} +struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, IC} + initial_condition :: IC mass :: ARRAY1D # [particle] radius :: ARRAY1D # [particle] elastic_modulus :: ELTYPE @@ -53,7 +53,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ST} <: SolidSystem{NDIMS} throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - return new{NDIMS, ELTYPE, typeof(mass), + return new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(source_terms)}(initial_condition, mass, radius, elastic_modulus, poissons_ratio, normal_stiffness, damping_coefficient, acceleration_, source_terms) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index fe949b9f1..5b3021b1c 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -49,7 +49,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method. where `beam` and `fixed_particles` are of type `InitialCondition`. """ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D, - K, PF, ST} <: SolidSystem{NDIMS} + K, PF, ST} <: SolidSystem{NDIMS, IC} initial_condition :: IC initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle] current_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle] diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index 3aa0af76d..70b179e3b 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -149,27 +149,35 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} particle_spacing, n_particles_per_dim) - if state_equation !== nothing - # Use hydrostatic pressure gradient and calculate density from inverse state - # equation, so don't pass fluid density. - fluid = RectangularShape(particle_spacing, n_particles_per_dim, zeros(NDIMS); - velocity, pressure, acceleration, state_equation, - mass=fluid_mass) - else - fluid = RectangularShape(particle_spacing, n_particles_per_dim, zeros(NDIMS); - density=fluid_density, velocity, pressure, - acceleration, state_equation, mass=fluid_mass) - end - boundary = InitialCondition(coordinates=boundary_coordinates, velocity=boundary_velocities, mass=boundary_masses, density=boundary_densities, particle_spacing=boundary_spacing) # Move the tank corner in the negative coordinate directions to the desired position - fluid.coordinates .+= min_coordinates boundary.coordinates .+= min_coordinates + if norm(fluid_size) > eps() + if state_equation !== nothing + # Use hydrostatic pressure gradient and calculate density from inverse state + # equation, so don't pass fluid density. + fluid = RectangularShape(particle_spacing, n_particles_per_dim, + zeros(NDIMS); + velocity, pressure, acceleration, state_equation, + mass=fluid_mass) + else + fluid = RectangularShape(particle_spacing, n_particles_per_dim, + zeros(NDIMS); + density=fluid_density, velocity, pressure, + acceleration, state_equation, mass=fluid_mass) + end + # Move the tank corner in the negative coordinate directions to the desired position + fluid.coordinates .+= min_coordinates + else + # Fluid is empty + fluid = InitialCondition(coordinates=zeros(ELTYPE, NDIMS, 0), density=1.0) + end + return new{NDIMS, 2 * NDIMS, ELTYPE}(fluid, boundary, fluid_size_, tank_size_, faces, face_indices, particle_spacing, spacing_ratio, n_layers, diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 6a41db723..0f0032607 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -114,6 +114,60 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/dam_break_2d_surface_tension.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d_surface_tension.jl"), + tspan=(0.0, 0.1)) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "sphere_surface_tension_2d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/sphere_surface_tension_3d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "sphere_surface_tension_3d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/falling_water_spheres_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "falling_water_spheres_2d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/falling_water_spheres_3d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "falling_water_spheres_3d.jl")) [ + r"┌ Info: The desired tank length in x-direction .*\n", + r"└ New tank length in x-direction.*\n", + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n", + r"┌ Info: The desired tank length in z-direction .*\n", + r"└ New tank length in z-direction.*\n", + ] + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/sphere_surface_tension_wall_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "sphere_surface_tension_wall_2d.jl")) + end + @trixi_testset "fluid/moving_wall_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), joinpath(examples_dir(), "fluid", @@ -190,7 +244,9 @@ @trixi_testset "n_body/n_body_benchmark_trixi.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "n_body", - "n_body_benchmark_trixi.jl")) + "n_body_benchmark_trixi.jl")) [ + r"WARNING: Method definition interact!.*\n", + ] end @trixi_testset "n_body/n_body_benchmark_reference.jl" begin diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 8a74de0ab..633a0e240 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -1,8 +1,8 @@ # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "Semidiscretization" begin # Mock systems - struct System1 <: TrixiParticles.System{3} end - struct System2 <: TrixiParticles.System{3} end + struct System1 <: TrixiParticles.System{3, Nothing} end + struct System2 <: TrixiParticles.System{3, Nothing} end system1 = System1() system2 = System2() @@ -39,7 +39,7 @@ struct BoundaryModelMock end # Mock fluid system - struct FluidSystemMock <: TrixiParticles.FluidSystem{2} end + struct FluidSystemMock <: TrixiParticles.FluidSystem{2, Nothing} end kernel = Val(:smoothing_kernel) Base.ndims(::Val{:smoothing_kernel}) = 2 diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index bf7c3720c..105fb0ee8 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -1,3 +1,4 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") +include("surface_tension.jl") diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl new file mode 100644 index 000000000..9758463f5 --- /dev/null +++ b/test/schemes/fluid/surface_tension.jl @@ -0,0 +1,92 @@ + +@testset verbose=true "Surface Tension" begin + @testset verbose=true "`cohesion_force_akinci`" begin + surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) + support_radius = 1.0 + m_b = 1.0 + pos_diff = [1.0, 1.0] + + # These values can be extracted from the graphs in the paper by Akinci et al. or by manual calculation. + # Additional digits have been accepted from the actual calculation. + test_distance = 0.1 + val = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance) * test_distance + @test isapprox(val[1], 0.1443038770421044, atol=6e-15) + @test isapprox(val[2], 0.1443038770421044, atol=6e-15) + + # Maximum repulsion force + test_distance = 0.01 + max = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance) * test_distance + @test isapprox(max[1], 0.15913517632298307, atol=6e-15) + @test isapprox(max[2], 0.15913517632298307, atol=6e-15) + + # Near 0 + test_distance = 0.2725 + zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance) * test_distance + @test isapprox(zero[1], 0.0004360543645195717, atol=6e-15) + @test isapprox(zero[2], 0.0004360543645195717, atol=6e-15) + + # Maximum attraction force + test_distance = 0.5 + maxa = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance) * test_distance + @test isapprox(maxa[1], -0.15915494309189535, atol=6e-15) + @test isapprox(maxa[2], -0.15915494309189535, atol=6e-15) + + # Should be 0 + test_distance = 1.0 + zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance) * test_distance + @test isapprox(zero[1], 0.0, atol=6e-15) + @test isapprox(zero[2], 0.0, atol=6e-15) + end + + @testset verbose=true "adhesion_force_akinci" begin + surface_tension = TrixiParticles.SurfaceTensionAkinci(surface_tension_coefficient=1.0) + support_radius = 1.0 + m_b = 1.0 + pos_diff = [1.0, 1.0] + + # These values can be extracted from the graphs in the paper by Akinci et al. or by manual calculation. + # Additional digits have been accepted from the actual calculation. + test_distance = 0.1 + zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance, 1.0) * + test_distance + @test isapprox(zero[1], 0.0, atol=6e-15) + @test isapprox(zero[2], 0.0, atol=6e-15) + + test_distance = 0.5 + zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance, 1.0) * + test_distance + @test isapprox(zero[1], 0.0, atol=6e-15) + @test isapprox(zero[2], 0.0, atol=6e-15) + + # Near 0 + test_distance = 0.51 + zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance, 1.0) * + test_distance + @test isapprox(zero[1], -0.002619160170741761, atol=6e-15) + @test isapprox(zero[2], -0.002619160170741761, atol=6e-15) + + # Maximum adhesion force + test_distance = 0.75 + max = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance, 1.0) * + test_distance + @test isapprox(max[1], -0.004949747468305833, atol=6e-15) + @test isapprox(max[2], -0.004949747468305833, atol=6e-15) + + # Should be 0 + test_distance = 1.0 + zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, test_distance, 1.0) * + test_distance + @test isapprox(zero[1], 0.0, atol=6e-15) + @test isapprox(zero[2], 0.0, atol=6e-15) + end +end diff --git a/test/systems/boundary_system.jl b/test/systems/boundary_system.jl index dad30e71c..2ffa40af4 100644 --- a/test/systems/boundary_system.jl +++ b/test/systems/boundary_system.jl @@ -111,7 +111,7 @@ system = BoundarySPHSystem(initial_condition, model) - show_compact = "BoundarySPHSystem{2}((hydrodynamic_mass = 3,), nothing) with 2 particles" + show_compact = "BoundarySPHSystem{2}((hydrodynamic_mass = 3,), nothing, 0.0) with 2 particles" @test repr(system) == show_compact show_box = """ @@ -121,6 +121,7 @@ │ #particles: ………………………………………………… 2 │ │ boundary model: ……………………………………… (hydrodynamic_mass = 3,) │ │ movement function: ……………………………… nothing │ + │ adhesion coefficient: ……………………… 0.0 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 0749261bd..68f69aea8 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -193,7 +193,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -206,6 +206,7 @@ │ smoothing kernel: ………………………………… Val │ │ viscosity: …………………………………………………… nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ + │ surface tension: …………………………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ source terms: …………………………………………… Nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""