svchb · svchb · Jul 31, 2026
diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md
@@ -220,9 +220,9 @@ Pages = [joinpath("general", "corrections.jl")]
 
 ### Overview of surface normal calculation in SPH
 
-Surface normals are essential for modeling surface tension as they provide the directionality
-of forces acting at the fluid interface. They are calculated based on the particle properties and
-their spatial distribution.
+Surface normals provide the directionality of forces acting at the fluid interface. They are
+used by the full Akinci model and both Morris models, but not by the cohesion-only Akinci model.
+They are calculated based on the particle properties and their spatial distribution.
 
 #### Color field and gradient-based surface normals
 
@@ -296,6 +296,32 @@ In the following table some values are shown for reference. The values marked wi
 | **Water**       | 0.07288  [Lange](@cite Lange2005)               |
 | **Mercury**     | 0.486502 [Lange](@cite Lange2005)               |
 
+### Model configuration
+
+All surface tension coefficients must be finite and non-negative. A zero coefficient disables
+the fluid-fluid surface force. Wall adhesion is controlled independently by the boundary's
+`adhesion_coefficient`.
+
+`CohesionForceAkinci` only evaluates the pairwise cohesion and optional wall-adhesion forces.
+It does not require surface normals or `reference_particle_spacing`. The full
+`SurfaceTensionAkinci` model and both Morris models require a surface-normal method. When one
+of these models is selected without an explicit method, `ColorfieldSurfaceNormal()` is used.
+Runnable configurations are available in `examples/fluid/cohesion_force_akinci_2d.jl` and
+`examples/fluid/akinci_wetting_2d.jl`.
+
+!!! warning "Akinci coefficients in two dimensions"
+    The Akinci cohesion and adhesion kernels implemented here use the normalization published
+    for the three-dimensional model. In two-dimensional simulations, their coefficients are
+    empirical numerical parameters rather than resolution-independent physical values in N/m.
+    Recheck the coefficient when changing particle spacing or smoothing length.
+
+    Both published kernels have dimensions of inverse volume: the cohesion numerator scales
+    with ``L^6`` and its denominator with ``L^9``, while the fourth-root term in the adhesion
+    kernel changes its ``L^{-3.25}`` prefactor to ``L^{-3}``. Particle mass scales with
+    ``L^D`` in ``D`` dimensions. At a fixed smoothing-length-to-spacing ratio, the resulting
+    pair contribution therefore scales with ``L^{D-3}``: it is resolution-independent in 3D
+    but retains an inverse-length factor in 2D.
+
 ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf)
 
 The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components,

diff --git a/examples/fluid/akinci_wetting_2d.jl b/examples/fluid/akinci_wetting_2d.jl
@@ -0,0 +1,67 @@
+# ==========================================================================================
+# 2D Akinci Surface Tension and Wall Adhesion
+#
+# A circular drop rests on the bottom wall. Set `wetting=true` to compare stronger wall
+# adhesion with a surface-tension-dominated non-wetting setup. The full Akinci model uses
+# surface normals, while its wall interaction uses the Akinci adhesion kernel.
+# ==========================================================================================
+
+using TrixiParticles
+using OrdinaryDiffEqLowStorageRK
+
+particle_spacing = 0.005
+fluid_density = 1000.0
+sound_speed = 100.0
+gravity = 9.81
+tspan = (0.0, 0.3)
+wetting = false
+
+if wetting
+    surface_tension_coefficient = 0.01
+    adhesion_coefficient = 1.0
+    nu = 0.0005
+else
+    surface_tension_coefficient = 2.0
+    adhesion_coefficient = 0.001
+    nu = 0.001
+end
+
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=1)
+
+tank = RectangularTank(particle_spacing, (0.0, 0.0), (0.5, 0.1), fluid_density;
+                       n_layers=3, faces=(true, true, true, false),
+                       acceleration=(0.0, -gravity), state_equation)
+drop = SphereShape(particle_spacing, 0.05, (0.25, 0.05), fluid_density;
+                   sphere_type=VoxelSphere())
+
+smoothing_length = particle_spacing - eps()
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+alpha = 8 * nu / (smoothing_length * sound_speed)
+viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
+
+fluid_system = WeaklyCompressibleSPHSystem(drop; smoothing_kernel, smoothing_length,
+                                           density_calculator=ContinuityDensity(),
+                                           state_equation, viscosity,
+                                           acceleration=(0.0, -gravity),
+                                           surface_tension=SurfaceTensionAkinci(;
+                                                                                surface_tension_coefficient),
+                                           correction=AkinciFreeSurfaceCorrection(fluid_density),
+                                           reference_particle_spacing=particle_spacing)
+
+boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system,
+                                             boundary_density_calculator=AdamiPressureExtrapolation(),
+                                             viscosity=ViscosityAdami(; nu=4 * nu),
+                                             clip_negative_pressure=true)
+boundary_system = WallBoundarySystem(tank.boundary, boundary_model;
+                                     adhesion_coefficient)
+
+semi = Semidiscretization(fluid_system, boundary_system)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+saving_callback = SolutionSavingCallback(dt=0.01)
+callbacks = CallbackSet(info_callback, saving_callback)
+
+sol = solve(ode, RDPK3SpFSAL35(); abstol=1e-7, reltol=1e-4,
+            save_everystep=false, callback=callbacks)
diff --git a/examples/fluid/cohesion_force_akinci_2d.jl b/examples/fluid/cohesion_force_akinci_2d.jl
@@ -0,0 +1,56 @@
+# ==========================================================================================
+# 2D Cohesion-Only Akinci Surface Force
+#
+# This example evolves a rectangular fluid patch with `CohesionForceAkinci`. The
+# cohesion-only model does not calculate surface normals and therefore does not require
+# `reference_particle_spacing`. In 2D, its coefficient is an empirical numerical parameter.
+# ==========================================================================================
+
+using TrixiParticles
+using OrdinaryDiffEqLowStorageRK
+
+particle_spacing = 0.025
+fluid_size = (0.2, 0.1)
+fluid_density = 1000.0
+sound_speed = 20.0
+tspan = (0.0, 0.2)
+
+fluid = RectangularShape(particle_spacing,
+                         round.(Int, fluid_size ./ particle_spacing),
+                         zeros(length(fluid_size)); density=fluid_density)
+
+smoothing_length = particle_spacing - eps()
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7, clip_negative_pressure=true)
+
+nu = 0.01
+alpha = 8 * nu / (smoothing_length * sound_speed)
+viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
+surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.001)
+
+fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length,
+                                           density_calculator=SummationDensity(),
+                                           state_equation, viscosity, surface_tension,
+                                           source_terms=SourceTermDamping(damping_coefficient=0.5))
+
+semi = Semidiscretization(fluid_system)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+saving_callback = SolutionSavingCallback(dt=0.02)
+stepsize_callback = StepsizeCallback(cfl=0.5)
+callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, save_everystep=false, callback=callbacks)
+
+v_ode, u_ode = sol.u[end].x
+v = TrixiParticles.wrap_v(v_ode, fluid_system, semi)
+velocity = TrixiParticles.current_velocity(v, fluid_system)
+total_momentum = vec(sum(velocity .* transpose(fluid_system.mass); dims=2))
+center_of_mass_velocity = total_momentum / sum(fluid_system.mass)
+final_kinetic_energy = kinetic_energy(fluid_system, nothing, nothing, v_ode, u_ode, semi,
+                                      sol.t[end])
+
+@info "Cohesion diagnostics" center_of_mass_velocity final_kinetic_energy
diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -65,12 +65,15 @@ viscosity_fluid = ArtificialViscosityMonaghan(; alpha, beta=0.0)
 density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 # density_diffusion = DensityDiffusionAntuono(delta=0.1)
 
+surface_tension = nothing
+reference_particle_spacing = 0
+
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothing_length,
                                            density_calculator=fluid_density_calculator,
                                            state_equation, viscosity=viscosity_fluid,
                                            density_diffusion, acceleration=(0.0, -gravity),
-                                           correction=nothing, surface_tension=nothing,
-                                           reference_particle_spacing=0)
+                                           correction=nothing, surface_tension,
+                                           reference_particle_spacing)
 
 # ==========================================================================================
 # ==== Boundary

diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl
@@ -33,12 +33,15 @@ nu_sim_water = nu_ratio * nu_sim_oil
 
 oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 
-# TODO: broken if both systems use surface tension
+# A physically consistent two-phase surface-tension interaction requires an explicit
+# interface model. Until that model is available, this example focuses on density and
+# viscosity contrast and keeps surface tension disabled on both fluids.
+surface_tension = nothing
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl");
               sol=nothing, fluid_particle_spacing, tspan,
               viscosity_fluid=ViscosityMorris(nu=nu_sim_water), smoothing_length, gravity,
               density_diffusion=nothing, sound_speed, prefix="",
-              reference_particle_spacing=fluid_particle_spacing)
+              surface_tension)
 
 # ==========================================================================================
 # ==== Setup oil layer
@@ -65,18 +68,7 @@ oil_system = WeaklyCompressibleSPHSystem(oil;
                                          state_equation=oil_eos,
                                          viscosity=oil_viscosity,
                                          acceleration=(0.0, -gravity),
-                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01),
-                                         correction=AkinciFreeSurfaceCorrection(oil_density),
-                                         reference_particle_spacing=fluid_particle_spacing)
-
-# oil_system = WeaklyCompressibleSPHSystem(oil;
-#                                          smoothing_kernel, smoothing_length,
-#                                          density_calculator=fluid_density_calculator,
-#                                          state_equation=oil_eos,
-#                                          viscosity=oil_viscosity,
-#                                          acceleration=(0.0, -gravity),
-#                                          surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03),
-#                                          reference_particle_spacing=fluid_particle_spacing)
+                                         surface_tension)
 
 # ==========================================================================================
 # ==== Simulation

diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl
@@ -1,9 +1,8 @@
 # ==========================================================================================
-# 2D Falling Water Spheres Simulation (With and Without Surface Tension)
+# 2D Falling Water Spheres Simulation with Surface Tension
 #
 # This example simulates two circular water "spheres" falling under gravity.
-# One sphere includes a surface tension model (Akinci et al.), while the other does not.
-# This demonstrates the effect of surface tension on fluid behavior.
+# Both spheres belong to one fluid system and use the same surface tension model.
 # ==========================================================================================
 
 using TrixiParticles
@@ -44,6 +43,7 @@ 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))
+falling_spheres = isnothing(sphere2) ? sphere1 : union(sphere1, sphere2)
 
 # ==========================================================================================
 # ==== Fluid
@@ -53,28 +53,20 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,
 fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
 fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
-fluid_density_calculator = ContinuityDensity()
-
 nu = 0.005
 alpha = 8 * nu / (fluid_smoothing_length * sound_speed)
 viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
-density_diffusion = DensityDiffusionAntuono(delta=0.1)
 surface_tension_coefficient = 0.05
 surface_tension = SurfaceTensionAkinci(; surface_tension_coefficient)
+reference_particle_spacing = isnothing(surface_tension) ? 0 : fluid_particle_spacing
 
-sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1;
+sphere_surface_tension = EntropicallyDampedSPHSystem(falling_spheres;
                                                      smoothing_kernel=fluid_smoothing_kernel,
                                                      smoothing_length=fluid_smoothing_length,
                                                      sound_speed, viscosity,
                                                      density_calculator=ContinuityDensity(),
                                                      acceleration, surface_tension,
-                                                     reference_particle_spacing=fluid_particle_spacing)
-
-sphere = WeaklyCompressibleSPHSystem(sphere2; smoothing_kernel=fluid_smoothing_kernel,
-                                     smoothing_length=fluid_smoothing_length,
-                                     density_calculator=fluid_density_calculator,
-                                     state_equation, viscosity, density_diffusion,
-                                     acceleration)
+                                                     reference_particle_spacing)
 
 # ==========================================================================================
 # ==== Boundary
@@ -94,7 +86,7 @@ boundary_system = WallBoundarySystem(tank.boundary, boundary_model;
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system)
+semi = Semidiscretization(sphere_surface_tension, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=1000)

diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl
@@ -1,10 +1,8 @@
 # ==========================================================================================
-# 3D Falling Water Spheres Simulation (With and Without Surface Tension)
+# 3D Falling Water Spheres Simulation with Surface Tension
 #
 # This example extends `falling_water_spheres_2d.jl` to three dimensions.
-# It simulates two spherical volumes of water falling under gravity.
-# One sphere includes a surface tension model, while the other does not,
-# demonstrating the effect of surface tension in 3D.
+# It simulates two spherical volumes of water in one fluid system falling under gravity.
 # ==========================================================================================
 
 using TrixiParticles

diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl
@@ -23,7 +23,7 @@ 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`
+# The surface tension configurations below use a compact support of `2 * particle_spacing`.
 # smoothing_length = particle_spacing
 # smoothing_kernel = WendlandC2Kernel{2}()
 # nu = 0.01

diff --git a/examples/fluid/sphere_surface_tension_3d.jl b/examples/fluid/sphere_surface_tension_3d.jl
@@ -17,7 +17,7 @@ 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`
+# The surface tension configurations below use a compact support of `2 * particle_spacing`.
 smoothing_length = 1.0 * particle_spacing
 
 nu = 0.04

diff --git a/examples/fluid/sphere_surface_tension_wall_2d.jl b/examples/fluid/sphere_surface_tension_wall_2d.jl
@@ -63,6 +63,6 @@ sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1;
                                                      reference_particle_spacing=fluid_particle_spacing)
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl");
-              sphere=nothing, sphere1, adhesion_coefficient=0.001, wall_viscosity=4.0 * nu,
+              sphere1, sphere2=nothing, adhesion_coefficient=0.001, wall_viscosity=4.0 * nu,
               alpha, sound_speed, fluid_density, nu, fluid_particle_spacing, tspan,
               tank_size, fluid_smoothing_length, sphere_surface_tension)
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -33,7 +33,7 @@ Boundary model for [`WallBoundarySystem`](@ref).
                                 in areas of low pressure, against which the particle
                                 shifting technique is fighting.
 - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
+                                which is needed when using a surface-normal method.
 # Examples
 ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1)
 # Free-slip condition