Implement Morris surface tension model

NEWS.md

@@ -3,6 +3,18 @@
 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.7
+
+### Features
+
+- Adds the classic **Continuum Surface Force (CSF)** model based on Morris 2000 (#584), which computes
+  surface tension as a **body force** proportional to curvature and directed along the interface normal.
+  This method is efficient and accurate for capillary effects but does not explicitly conserve momentum.
+
+- Added the classic **Continuum Surface Stress (CSS)** model based on Morris 2000 (#584), which is
+  a **momentum-conserving** approach that formulates surface tension as the **divergence of a stress tensor**.
+  However, it requires additional computation and stabilization to handle **high-density interfaces** and reduce numerical instabilities.
+
 ## Version 0.2.6
 
 ### Features

Project.toml

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

docs/src/systems/fluid.md

@@ -1,11 +1,51 @@
 
 # [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).  
+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.
+Viscosity is a critical physical property governing momentum diffusion within a fluid.
+In the context of SPH, viscosity determines how rapidly velocity gradients are smoothed out,
+influencing key flow characteristics such as boundary layer formation, vorticity diffusion,
+and dissipation of kinetic energy. It also helps determine whether a flow is laminar or turbulent
+under a given set of conditions.
+
+Implementing viscosity correctly in SPH is essential for producing physically accurate results,
+and different methods exist to capture both numerical stabilization and true viscous effects.
+
+### Artificial (numerical) viscosity
+
+- Goal: Stabilize the simulation, capture shocks, and prevent unphysical particle interpenetration.
+- Method: Adds a dissipative (artificial) term to the momentum equations.
+- Typical Use: High-speed flows with strong shocks, astrophysical simulations, or situations where numerical damping is needed for stability.
+
+### Physical (real) viscosity
+
+- Goal: Model the actual viscous stresses of a fluid, aligned with a target Reynolds number or experimentally measured fluid properties.
+- Method: Introduces a force consistent with the Navier–Stokes viscous stress term.
+- Typical Use: Low-speed, incompressible or weakly compressible flows where matching real fluid behavior is important.
+
+### Model comparison
+
+#### ArtificialViscosityMonaghan
+
+- Best For: Compressible/high-speed flows, shock capturing, general purpose damping.
+- If you need: Stability in challenging flow regimes with potentially large density/pressure variations.
+
+#### ViscosityMorris
+
+- Best For: Moderate to low Mach number flows where realistic viscous behavior is desired.
+- If you need: Straightforward approach to physical viscosity that still works well in weakly compressible scenarios.
+
+#### ViscosityAdami
+
+- Best For: Incompressible or weakly compressible flows requiring accurate shear stress treatment.
+- If you need: Good boundary layer representation and accurate laminar flow with minimal compressibility effects.
+
+
+### API
 
 ```@autodocs
 Modules = [TrixiParticles]
@@ -19,40 +59,44 @@ Modules = [TrixiParticles]
 Pages = [joinpath("general", "corrections.jl")]
 ```
 
-
 ---
 
 ## [Surface Normals](@id surface_normal)
 
-### Overview of Surface Normal Calculation in SPH
+### 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.
+Surface normals are essential for modeling surface tension as they provide the directionalit
+of forces acting at the fluid interface. They are calculated based on the particle properties and
+their spatial distribution.
 
-#### Color Field and Gradient-Based Surface Normals
+#### 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 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:
+The simplest SPH formulation for a 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 \).
+- ``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
+#### 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
+#### 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:
+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.
 
@@ -65,17 +109,38 @@ 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.
+
+The surface tension coefficient ``\sigma`` is a physical parameter that quantifies the energy required to increase the surface area
+of a fluid by a unit amount. A higher value of ``\sigma`` indicates that the fluid resists changes to its surface area more strongly,
+causing droplets or bubbles to assume shapes (often spherical) that minimize their surface. In practice, ``\sigma`` can be measured
+experimentally through techniques such as the pendant drop method, the Wilhelmy plate method, or the du Noüy ring method,
+each of which relates a measurable force or change in shape to the fluid’s surface tension. For pure substances,
+tabulated reference values of ``\sigma`` at given temperatures are commonly used, while for mixtures or complex fluids,
+direct experimental measurements or specialized equations of state may be necessary.
+
+| **Fluid**    | **Surface Tension (``\sigma``) [N/m at 20°C]** |
+|--------------|----------------------------------------------:|
+| **Water**    | 0.0728                                        |
+| **Mercury**  | 0.485                                         |
+| **Ethanol**  | 0.0221                                        |
+| **Acetone**  | 0.0237                                        |
+| **Glycerol** | 0.0634                                        |
+| **Olive Oil**| ~0.032                                        |
+| **Gasoline** | ~0.022                                        |
+| **Mineral Oil** | ~0.030                                     |
 
-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.
+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.
@@ -97,15 +162,17 @@ C(r)=\frac{32}{\pi h_c^9}
 
 #### 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:
+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:
+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},
 ```
@@ -120,32 +187,57 @@ A(r) = \frac{0.007}{h_c^{3.25}}
 
 ---
 
+### Morris Surface Tension Model
+
+The method estimates curvature by combining particle color gradients and smoothing functions to derive surface normals.
+The computed curvature is then used to determine forces acting perpendicular to the interface.
+While this method provides accurate surface tension forces, it does not explicitly conserve momentum.
+
+In the Morris model, surface tension is computed based on local interface curvature ``\kappa`` and the unit surface normal ``\hat{n}``.
+By estimating ``\hat{n}`` and ``\kappa`` at each particle near the interface, the surface tension force for particle a can be written as:
+
+```math
+F_{\text{surface tension}} = - \sigma \frac{\kappa_a}{\rho_a}\hat{n}_a
+```
+This formulation focuses directly on geometric properties of the interface, making it relatively straightforward to implement when a reliable interface detection 
+(e.g., a color function) is available. However, accurately estimating ``\kappa`` and ``n`` may require fine resolutions.
+---
+
 ### 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.
+In addition to the simpler curvature-based formulation, Morris (2000) introduced a momentum-conserving approach.
+This method treats surface tension forces as arising from the divergence of a stress tensor, ensuring exact conservation
+of linear momentum and offering more robust behavior for high-resolution or long-duration simulations
+where accumulated numerical error can be significant.
 
 #### Stress Tensor Formulation
 
-The force is calculated as:
+The surface tension force can be seen as a divergence of a stress tensor ``S``
 ```math
 F_{\text{surface tension}} = \nabla \cdot S,
 ```
-where \( S \) is the stress tensor:
+with ``S`` defined as
 ```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.
+- ``\delta_s``: Surface delta function,
+- ``\hat{n}``: Unit normal vector,
+- ``I``: Identity matrix.
+
+That can be calculated as
+```math 
+\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab} 
+```
 
 #### 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.
+While momentum conservation makes this model attractive, it requires additional computational effort and stabilization
+techniques to address instabilities in high-density regions.
+
+### API
 
 ```@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.
+```

examples/fluid/falling_water_spheres_3d.jl

@@ -34,4 +34,4 @@ trixi_include(@__MODULE__,
               fluid_smoothing_length=fluid_smoothing_length,
               fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(),
               nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed),
-              surface_tension_coefficient=10, adhesion_coefficient=0.1)
+              surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05))

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -48,6 +48,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used)
 - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
                                 which currently is only needed when using surface tension.
+- `color_value`:                The value used to initialize the color of particles in the system.
 
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,