Implement Morris surface tension model

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

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"

docs/make.jl

@@ -140,8 +140,8 @@ makedocs(sitename="TrixiParticles.jl",
                      "Util" => joinpath("general", "util.md")
                  ],
                  "Systems" => [
-                     "Discrete Element Method (Solid)" => joinpath("systems",
-                                                                   "dem.md"),
+                     "Fluid Models" => joinpath("systems", "fluid.md"),
+                     "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",

docs/src/refs.bib

@@ -491,6 +491,19 @@ @Article{Morris1997
   publisher = {Elsevier BV},
 }
 
+@article{Morris2000,
+  author = {Morris, Joseph P.},
+  title = {Simulating surface tension with smoothed particle hydrodynamics},
+  journal = {International Journal for Numerical Methods in Fluids},
+  volume = {33},
+  number = {3},
+  pages = {333-353},
+  keywords = {interfacial flow, meshless methods, surface tension},
+  doi = {https://doi.org/10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7},
+  year = {2000}
+}
+
+
 @InProceedings{Mueller2003,
   author    = {M{\"u}ller, Matthias and Charypar, David and Gross, Markus},
   booktitle = {Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation},

docs/src/systems/fluid.md

@@ -0,0 +1,151 @@
+
+# [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)
+
+### 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 within the smoothed particle hydrodynamics (SPH) framework.
+
+#### Color Field and Gradient-Based Surface Normals
+
+The surface normal at a particle is derived from the color field, a scalar field assigned to particles to distinguish between different fluid phases or between fluid and air. The color field gradients point towards the interface, and the normalized gradient defines the surface normal direction.
+
+The simplest SPH formulation for surface normal, \( n_a \), is given as:
+```math
+n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab},
+```
+where:
+- \( c_b \) is the color field value for particle \( b \),
+- \( m_b \) is the mass of particle \( b \),
+- \( \rho_b \) is the density of particle \( b \),
+- \( \nabla_a W_{ab} \) is the gradient of the smoothing kernel \( W_{ab} \) with respect to particle \( a \).
+
+#### Normalization of Surface Normals
+
+The calculated normals are normalized to unit vectors:
+```math
+\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}.
+```
+Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.
+
+#### Handling Noise and Errors in Normal Calculation
+
+In regions distant from the interface, the calculated normals may be small or inaccurate due to the smoothing kernel's support radius. To mitigate this:
+1. Normals below a threshold are excluded from further calculations.
+2. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "surface_normal_sph.jl")]
+```
+
+---
+
+## [Surface Tension](@id surface_tension)
+
+### Introduction to Surface Tension in SPH
+
+Surface tension is a key phenomenon in fluid dynamics, influencing the behavior of droplets, bubbles, and fluid interfaces. In SPH, surface tension is modeled as forces arising due to surface curvature and particle interactions, ensuring realistic simulation of capillary effects, droplet coalescence, and fragmentation.
+
+### Akinci-Based Intra-Particle Force Surface Tension and Wall Adhesion Model
+
+The Akinci model divides surface tension into distinct force components:
+
+#### Cohesion Force
+
+The cohesion force captures the attraction between particles at the fluid interface, creating the effect of surface tension. It is defined by the distance between particles and the support radius \( h_c \), using a kernel-based formulation.
+
+**Key Features:**
+- Particles within half the support radius experience a repulsive force to prevent clustering.
+- Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion.
+
+Mathematically:
+```math
+F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},
+```
+where \( C(r) \), the cohesion kernel, 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
+
+The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. It acts based on the difference in surface normals:
+```math
+F_{\text{curvature}} = -\sigma (n_a - n_b),
+```
+where \( n_a \) and \( n_b \) are the surface normals of the interacting particles.
+
+#### Wall Adhesion Force
+
+This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. It uses a custom kernel with a peak at 0.75 times the support radius:
+```math
+F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},
+```
+where \( A(r) \) is the adhesion kernel:
+```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}
+```
+
+---
+
+### Morris-Based Momentum-Conserving Surface Tension Model
+
+In addition to the Akinci model, Morris (2000) introduced a momentum-conserving approach to surface tension. This model uses stress tensors to ensure exact conservation of linear momentum, providing a robust method for high-resolution simulations.
+
+#### Stress Tensor Formulation
+
+The force is calculated as:
+```math
+F_{\text{surface tension}} = \nabla \cdot S,
+```
+where \( S \) is the stress tensor:
+```math
+S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}),
+```
+with:
+- \( \delta_s \): Surface delta function,
+- \( \hat{n} \): Unit normal vector,
+- \( I \): Identity matrix.
+
+#### Advantages and Limitations
+
+While momentum conservation makes this model attractive, it requires additional computational effort and stabilization techniques to address instabilities in high-density regions.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "surface_tension.jl")]
+```
+
+This extended documentation provides a comprehensive view of the theoretical foundations and practical implementations of surface tension and surface normal calculations in SPH models.

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

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

examples/fluid/dam_break_oil_film_2d.jl

@@ -25,18 +25,31 @@ nu_ratio = nu_water / nu_oil
 nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
 nu_sim_water = nu_ratio * nu_sim_oil
 
+oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
+# Physical values
+nu_water = 8.9E-7
+nu_oil = 6E-5
+nu_ratio = nu_water / nu_oil
+
+nu_sim_oil = max(0.01 * smoothing_length * sound_speed, nu_oil)
+nu_sim_water = nu_ratio * nu_sim_oil
+
 oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 
+# TODO: broken if both are set to surface tension
 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, tspan=tspan, 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 +63,19 @@ 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)
+
+# oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
+#                                          oil_eos, smoothing_kernel,
+#                                          smoothing_length, viscosity=oil_viscosity,
+#                                          acceleration=(0.0, -gravity),
+#                                          surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03),
+#                                          reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Simulation