Generalize surface normal calc (#539)

* fix merge * fix merge * fix merge * fix merge * fix merge * more merge fixes * format * fix * fix * fix * fix * fix * more fixes * fix * fix * fix * more fixes * fix * fix * format * running * update * update * optimization * update * update * update * update * format * up * update * remove unused code * switch to support radius * update * fix doc * add adhesion coefficient * update_broken * update * format * add to examples * fix tests * fix merge * fix typos * add basic test * format * rename * format * rename container system * reduce example run time * update news and readme * format * fix test * reduce run time * format * fix * correct merge * update doc test * revert * fix tests * fix * fix * review comments * fix * fix * correct some stuff * init with empty initial condition * review update * rename function * remove unnecessary if * docs * fix doc * revert one change * fix typo * update * try to avoid allocs * fix mem allocs * review update * update docs * update * review * fix * format * make examples smaller * reduce resolution * review * rename * use trixi_include * format * use trixi_include * update * update * review update * format * fix test errors * update * fix * remove invalid surface normals * add * fix * update * update * update * fix * update * update * update * remove unused * remove unused * format * fix tests * fix * fix * fix naming * format * fix test * set test up for 1.11 * fix test * rename to number_density * back merge * fix * format * fix * fix tests * back merge from #584 * move additonal changes from #584 * add test * add test * format * Increase errors for 1.11 * Fix invalidations * Fix tests * Update ci.yml * revert * Update ci.yml * Update test/validation/validation.jl Co-authored-by: Erik Faulhaber <[email protected]> * format * Update news and set to 0.2.4 * update * fix MD format * Update NEWS.md * some fixes * add docs * review fixes * remove independent setting of smoothing_kernel and smoothing_length * remove calls to surface_normal method smoothing kernel * format * update news * typo * Update NEWS.md * forgot some renames * fix doc tests * fix docs * fix tests * format * Update NEWS.md * fix tests * Update NEWS.md * fix * fix * fix * fix * review fixes * review comments * fix test * typo * review updates * review update * add boundary system * update * update * format * remove color since it is unused in this PR * review comment * format * merge error * review comment * fix * remove surface_normal_method() * merge error * function was renamed in PN * format * fix test * fix error * fix * review comments * review comments * review comments * format * merge * format * add comment * remove ideal_neighbor() * format * typo * add missing docs * only initialize necessary values * fix test * fix examples --------- Co-authored-by: Erik Faulhaber <[email protected]>
trixi-framework · Jan 30, 2025 · 2540f08 · svchb · Jan 30, 2025 · 2540f08
1 parent 3d73412
commit 2540f08
Show file tree

Hide file tree

Showing 28 changed files with 788 additions and 273 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -3,9 +3,20 @@
 TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.2.6
+
+### Features
+
+- Support for surface tension was added to EDAC (#539)
+
+### Refactored
+
+- Surface normal calculation was moved from surface_tension.jl to surface_normal_sph.jl (#539)
+
 ## Version 0.2.5
 
 ### Features
+
 - Add particle packing for 2D (.asc) and 3D (.stl) geometries (#529)
 
 ### Compatibility Changes
@@ -26,7 +37,6 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - A new user tutorial was added (#514)
 - Several docs fixes (#663, #659, #637, #658, #664)
 
-
 ## Version 0.2.3
 
 ### Highlights

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

diff --git a/docs/make.jl b/docs/make.jl
@@ -140,12 +140,14 @@ makedocs(sitename="TrixiParticles.jl",
                      "Util" => joinpath("general", "util.md")
                  ],
                  "Systems" => [
-                     "Discrete Element Method (Solid)" => joinpath("systems",
-                                                                   "dem.md"),
-                     "Weakly Compressible SPH (Fluid)" => joinpath("systems",
-                                                                   "weakly_compressible_sph.md"),
-                     "Entropically Damped Artificial Compressibility for SPH (Fluid)" => joinpath("systems",
-                                                                                                  "entropically_damped_sph.md"),
+                     "Fluid Models" => [
+                         "Overview" => joinpath("systems", "fluid.md"),
+                         "Weakly Compressible SPH (Fluid)" => joinpath("systems",
+                                                                       "weakly_compressible_sph.md"),
+                         "Entropically Damped Artificial Compressibility for SPH (Fluid)" => joinpath("systems",
+                                                                                                      "entropically_damped_sph.md")
+                     ],
+                     "Discrete Element Method (Solid)" => joinpath("systems", "dem.md"),
                      "Total Lagrangian SPH (Elastic Structure)" => joinpath("systems",
                                                                             "total_lagrangian_sph.md"),
                      "Boundary" => joinpath("systems", "boundary.md")

diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md
@@ -0,0 +1,89 @@
+# [Fluid Models](@id fluid_models)
+Currently available fluid methods are the [weakly compressible SPH method](@ref wcsph) and the [entropically damped artificial compressibility for SPH](@ref edac).  
+This page lists models and techniques that apply to both of these methods.  
+
+## [Viscosity](@id viscosity_wcsph)
+
+TODO: Explain viscosity.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "viscosity.jl")]
+```
+
+## [Corrections](@id corrections)
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("general", "corrections.jl")]
+```
+
+
+
+## [Surface Normals](@id surface_normal)
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "surface_normal_sph.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 support radius ``h_c`` and the distance between particles.
+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 particles from 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/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md
@@ -44,15 +44,6 @@ Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "state_equations.jl")]
 ```
 
-## [Viscosity](@id viscosity_wcsph)
-
-TODO: Explain viscosity.
-
-```@autodocs
-Modules = [TrixiParticles]
-Pages = [joinpath("schemes", "fluid", "viscosity.jl")]
-```
-
 ## [Density Diffusion](@id density_diffusion)
 
 Density diffusion can be used with [`ContinuityDensity`](@ref) to remove the noise in the
@@ -118,72 +109,3 @@ term is very expensive and adds about 40--50% of computational cost.
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "weakly_compressible_sph", "density_diffusion.jl")]
 ```
-
-## [Corrections](@id corrections)
-
-```@autodocs
-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
@@ -61,7 +61,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            smoothing_length, viscosity=viscosity,
                                            density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity), correction=nothing,
-                                           surface_tension=nothing)
+                                           surface_tension=nothing,
+                                           reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Boundary
@@ -70,7 +71,8 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
                                              state_equation=state_equation,
                                              boundary_density_calculator,
                                              smoothing_kernel, smoothing_length,
-                                             correction=nothing)
+                                             correction=nothing,
+                                             reference_particle_spacing=fluid_particle_spacing)
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0)
 

diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl
@@ -30,13 +30,16 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
               viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
-              gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed)
+              gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed,
+              prefix="")
 
 # ==========================================================================================
 # ==== Setup oil layer
 
 oil_size = (W, 0.1 * H)
 oil_density = 700.0
+oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1,
+                            clip_negative_pressure=false)
 
 oil = RectangularShape(fluid_particle_spacing,
                        round.(Int, oil_size ./ fluid_particle_spacing),
@@ -50,12 +53,12 @@ end
 oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density,
                                        exponent=1, clip_negative_pressure=false)
 oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
-                                         oil_state_equation, smoothing_kernel,
+                                         oil_eos, smoothing_kernel,
                                          smoothing_length, viscosity=oil_viscosity,
-                                         density_diffusion=density_diffusion,
                                          acceleration=(0.0, -gravity),
-                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02),
-                                         correction=AkinciFreeSurfaceCorrection(oil_density))
+                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01),
+                                         correction=AkinciFreeSurfaceCorrection(oil_density),
+                                         reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Simulation

diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl
@@ -40,7 +40,10 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,
 
 # ==========================================================================================
 # ==== Fluid
-fluid_smoothing_length = 1.0 * fluid_particle_spacing
+
+# Using a smoothing_length of exactly 1.0 * fluid_particle is necessary for this model to be accurate.
+# This yields some numerical issues though which can be circumvented by subtracting eps().
+fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
 fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
@@ -50,13 +53,13 @@ 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,
+sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
                                                      fluid_smoothing_length,
-                                                     viscosity=viscosity,
+                                                     sound_speed, viscosity=viscosity,
+                                                     density_calculator=ContinuityDensity(),
                                                      acceleration=(0.0, -gravity),
                                                      surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
-                                                     correction=AkinciFreeSurfaceCorrection(fluid_density))
+                                                     reference_particle_spacing=fluid_particle_spacing)
 
 sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator,
                                      state_equation, fluid_smoothing_kernel,
@@ -72,19 +75,20 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
                                              state_equation=state_equation,
                                              boundary_density_calculator,
                                              fluid_smoothing_kernel, fluid_smoothing_length,
-                                             viscosity=ViscosityAdami(nu=wall_viscosity))
+                                             viscosity=ViscosityAdami(nu=wall_viscosity),
+                                             reference_particle_spacing=fluid_particle_spacing)
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
-                                    adhesion_coefficient=0.001)
+                                    adhesion_coefficient=1.0)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere)
+semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="",
-                                         write_meta_data=true)
+saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out",
+                                         prefix="", write_meta_data=true)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl
@@ -3,35 +3,35 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-fluid_particle_spacing = 0.008
+fluid_particle_spacing = 0.005
 
 # ==========================================================================================
 # ==== Experiment Setup
 gravity = 9.81
-nu = 0.01
+nu = 0.001
 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_center = (0.5, 0.5, 0.075)
+sphere2_center = (1.5, 0.5, 0.075)
 sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0))
 sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.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),
+              fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.1),
               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)
+              surface_tension_coefficient=10, adhesion_coefficient=0.1)