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 <44124897+efaulhaber@users.noreply.github.com>"]
-version = "0.2.7-dev"
+version = "0.2.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

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

@@ -1,10 +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]
@@ -18,72 +59,185 @@ 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 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
+
+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 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``.
+
+#### 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)
 
-### 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:
+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                                     |
+
+
+### 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:
-- ``\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
+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}
+(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
+#### 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)
+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
 
-#### 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
+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},
 ```
-with ``A`` being the adhesion kernel defined as
+where \( A(r) \) is the adhesion kernel:
 ```math
-A(r)= \frac{0.007}{h_c^{3.25}}
+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 \\
+\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 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 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 surface tension force can be seen as a divergence of a stress tensor ``S``
+```math
+F_{\text{surface tension}} = \nabla \cdot S,
+```
+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.
+
+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.
+
+### API
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "surface_tension.jl")]
-```
+```

examples/fluid/dam_break_oil_film_2d.jl

@@ -26,7 +26,17 @@ 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,
@@ -60,6 +70,13 @@ oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
                                          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
 semi = Semidiscretization(fluid_system, oil_system, boundary_system,

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