Merge branch 'main' of github.com:svchb/TrixiParticles.jlOpen

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.
 

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.1.2-dev"
+version = "0.1.3-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

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

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")]
+```

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

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))

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);

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)

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);

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)