Add simple SGS to Adami and Morris Viscosity

examples/fluid/dam_break_2d.jl

@@ -47,10 +47,21 @@ smoothing_length = 1.75 * fluid_particle_spacing
 smoothing_kernel = WendlandC2Kernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
-viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
-# nu = 0.02 * smoothing_length * sound_speed/8
-# viscosity = ViscosityMorris(nu=nu)
+alpha = 0.02
+viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
+# A typical formula to convert Artificial viscosity to a
+# kinematic viscosity is provided by Monaghan as
+# nu = alpha * smoothing_length * sound_speed/8
+
+# Alternatively a kinematic viscosity for water can be set
+# nu = 1.0e-6
+
+# This allows the use of a physical viscosity model like:
 # viscosity = ViscosityAdami(nu=nu)
+# or with additional dissipation through a Smagorinsky model
+# viscosity = ViscosityAdamiSGS(nu=nu)
+# For more details see the documentation [Viscosity model overview](@ref viscosity_sph).
+
 # Alternatively the density diffusion model by Molteni & Colagrossi can be used,
 # which will run faster.
 # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
@@ -67,12 +78,17 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
 # ==========================================================================================
 # ==== Boundary
 boundary_density_calculator = AdamiPressureExtrapolation()
+viscosity_wall = nothing
+# For a no-slip condition the corresponding wall viscosity without SGS can be set
+#viscosity_wall = ViscosityAdami(nu=nu)
+#viscosity_wall = ViscosityMorris(nu=nu)
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              boundary_density_calculator,
                                              smoothing_kernel, smoothing_length,
                                              correction=nothing,
-                                             reference_particle_spacing=0)
+                                             reference_particle_spacing=0,
+                                             viscosity=viscosity_wall)
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coefficient=0.0)
 

src/TrixiParticles.jl

@@ -73,7 +73,8 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel
 export StateEquationCole, StateEquationIdealGas
-export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris
+export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS,
+       ViscosityMorrisSGS
 export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
        DensityDiffusionAntuono
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,

src/io/write_vtk.jl

@@ -357,7 +357,9 @@ end
 
 write2vtk!(vtk, viscosity::Nothing) = vtk
 
-function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris})
+function write2vtk!(vtk,
+                    viscosity::Union{ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS,
+                                     ViscosityMorrisSGS})
     vtk["viscosity_nu"] = viscosity.nu
     vtk["viscosity_epsilon"] = viscosity.epsilon
 end

src/schemes/fluid/viscosity.jl

@@ -176,6 +176,32 @@ struct ViscosityAdami{ELTYPE}
     end
 end
 
+function adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
+                               m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+    eta_a = nu_a * rho_a
+    eta_b = nu_b * rho_b
+
+    eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b)
+
+    tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2)
+
+    volume_a = m_a / rho_a
+    volume_b = m_b / rho_b
+
+    # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001
+    # They argued that the formulation is more flexible because of the possibility to formulate
+    # different inter-particle averages or to assume different inter-particle distributions.
+    # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration.
+    #
+    # TODO: Is there a better formulation to discretize the Laplace operator?
+    # Because when using this formulation for the pressure acceleration, it is not
+    # energy conserving.
+    # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394
+    visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a
+
+    return visc .* v_diff
+end
+
 @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system,
                                              v_particle_system, v_neighbor_system,
                                              particle, neighbor, pos_diff,
@@ -185,6 +211,7 @@ end
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
     smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
                                viscosity_model(neighbor_system, particle_system),
@@ -197,36 +224,232 @@ end
     v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
     v_diff = v_a - v_b
 
-    eta_a = nu_a * rho_a
-    eta_b = nu_b * rho_b
+    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
+                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+end
 
-    eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b)
+function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length,
+                             sound_speed)
+    return viscosity.nu
+end
+
+@propagate_inbounds function viscous_velocity(v, system, particle)
+    return current_velocity(v, system, particle)
+end
+
+@doc raw"""
+    ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01)
+
+Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami)
+by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type closure.
+The effective kinematic viscosity is defined as
+
+```math
+\nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}},
+```
+
+with
+
+```math
+\nu_{\mathrm{SGS}} = (C_S * h)^2 * |S|,
+```
+
+and an approximation for the strain rate magnitude given by
+
+```math
+|S| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon},
+```
+
+where:
+- $C_S$ is the Smagorinsky constant (typically 0.1 to 0.2),
+- $h$ is the local smoothing length, and
+
+The effective dynamic viscosities are then computed as
+```math
+\eta_{a,\mathrm{eff}} = \rho_a\, \nu_{\mathrm{eff}},
+```
+and averaged as
+```math
+\bar{\eta}_{ab} = \frac{2 \eta_{a,\mathrm{eff}} \eta_{b,\mathrm{eff}}}{\eta_{a,\mathrm{eff}}+\eta_{b,\mathrm{eff}}}.
+```
 
+This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation.
+
+# Keywords
+- `nu`:      Standard kinematic viscosity.
+- `C_S`:     Smagorinsky constant.
+- `epsilon=0.01`: Parameter to prevent singularities
+"""
+struct ViscosityAdamiSGS{ELTYPE}
+    nu::ELTYPE      # kinematic viscosity [e.g., 1e-6 m²/s]
+    C_S::ELTYPE     # Smagorinsky constant [e.g., 0.1-0.2]
+    epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001]
+end
+
+# Convenient constructor with default values for C_S and epsilon.
+ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon)
+
+@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system,
+                                                            neighbor_system,
+                                                            v_particle_system,
+                                                            v_neighbor_system,
+                                                            particle, neighbor, pos_diff,
+                                                            distance, sound_speed, m_a, m_b,
+                                                            rho_a, rho_b, grad_kernel)
+    epsilon = viscosity.epsilon
+
+    smoothing_length_particle = smoothing_length(particle_system, particle)
+    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
     smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
-    tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2)
 
-    volume_a = m_a / rho_a
-    volume_b = m_b / rho_b
+    nu_a = kinematic_viscosity(particle_system,
+                               viscosity_model(neighbor_system, particle_system),
+                               smoothing_length_particle, sound_speed)
+    nu_b = kinematic_viscosity(neighbor_system,
+                               viscosity_model(particle_system, neighbor_system),
+                               smoothing_length_neighbor, sound_speed)
 
-    # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001
-    # They argued that the formulation is more flexible because of the possibility to formulate
-    # different inter-particle averages or to assume different inter-particle distributions.
-    # Ramachandran (2019) and Adami (2012) use this formulation also for the pressure acceleration.
+    v_a = viscous_velocity(v_particle_system, particle_system, particle)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_diff = v_a - v_b
+
+    # ------------------------------------------------------------------------------
+    # SGS part: Compute the subgrid-scale eddy viscosity.
+    # ------------------------------------------------------------------------------
+    # In classical LES the Smagorinsky model defines:
+    #   ν_SGS = (C_S Δ)^2 |S|,
+    # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ).
     #
-    # TODO: Is there a better formulation to discretize the Laplace operator?
-    # Because when using this formulation for the pressure acceleration, it is not
-    # energy conserving.
-    # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394
-    visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a
+    # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly.
+    # A common low-order surrogate is to approximate the strain‐rate magnitude by a
+    # finite-difference along each particle pair:
+    #
+    #   |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ),
+    #
+    # where δ regularizes the denominator to avoid singularities when particles are very close.
+    #
+    # This yields:
+    #   S_mag = norm(v_diff) / (distance + ε)
+    #
+    # and then the Smagorinsky eddy viscosity:
+    #   ν_SGS = (C_S * h̄)^2 * S_mag.
+    #
+    S_mag = norm(v_diff) / (distance + epsilon)
+    nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag
 
-    return visc .* v_diff
+    # Effective kinematic viscosity is the sum of the standard and SGS parts.
+    nu_a = nu_a + nu_SGS
+    nu_b = nu_b + nu_SGS
+
+    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
+                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length,
+function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length,
                              sound_speed)
     return viscosity.nu
 end
 
-@propagate_inbounds function viscous_velocity(v, system, particle)
-    return current_velocity(v, system, particle)
+@doc raw"""
+    ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001)
+
+Subgrid-scale (SGS) viscosity model based on the formulation by [Morris (1997)](@cite Morris1997),
+extended with a Smagorinsky-type eddy viscosity term for modeling turbulent flows.
+
+The acceleration on particle `a` due to viscosity interaction with particle `b` is calculated as:
+```math
+\frac{d v_a}{dt} = \sum_b m_b \frac{\mu_{a,\mathrm{eff}} + \mu_{b,\mathrm{eff}}}{\rho_a \rho_b} \frac{r_{ab} \cdot \nabla_a W_{ab}}{r_{ab}^2 + \epsilon h_{ab}^2} v_{ab}
+where $v_{ab} = v_a - v_b$, $r_{ab} = r_a - r_b$, $r_{ab} = \|r_{ab}\|$, $h_{ab} = (h_a + h_b)/2$, $W_{ab}$ is the smoothing kernel, $\rho$ is density, $m$ is mass, and $\epsilon$ is a regularization parameter.
+
+The effective dynamic viscosity $\mu_{i,\mathrm{eff}}$ for each particle `i` (a or b) includes both the standard molecular viscosity and an SGS eddy viscosity contribution:
+```math
+\mu_{i,\mathrm{eff}} = \rho_i \nu_{i,\mathrm{eff}} = \rho_i (\nu_{\mathrm{std}} + \nu_{i,\mathrm{SGS}})
+```
+The standard kinematic viscosity $\nu_{\mathrm{std}}$ is provided by the `nu` parameter. The SGS kinematic viscosity is calculated using the Smagorinsky model:
+```math
+\nu_{i,\mathrm{SGS}} = (C_S h_{ab})^2 |\bar{S}_{ab}|
+```
+where $C_S$ is the Smagorinsky constant. This implementation uses a simplified pairwise approximation for the strain rate tensor magnitude $|\bar{S}_{ab}|$:
+```math
+|\bar{S}_{ab}| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon}
+```
+
+This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation.
+
+# Keywords
+- `nu`:      Standard kinematic viscosity.
+- `C_S`:     Smagorinsky constant.
+- `epsilon=0.01`: Parameter to prevent singularities
+"""
+struct ViscosityMorrisSGS{ELTYPE}
+    nu::ELTYPE      # kinematic viscosity [e.g., 1e-6 m²/s]
+    C_S::ELTYPE     # Smagorinsky constant [e.g., 0.1-0.2]
+    epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001]
+end
+
+ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon)
+
+@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system,
+                                                             neighbor_system,
+                                                             v_particle_system,
+                                                             v_neighbor_system,
+                                                             particle, neighbor, pos_diff,
+                                                             distance, sound_speed, m_a,
+                                                             m_b, rho_a, rho_b, grad_kernel)
+    epsilon = viscosity.epsilon
+
+    smoothing_length_particle = smoothing_length(particle_system, particle)
+    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
+
+    nu_a = kinematic_viscosity(particle_system,
+                               viscosity_model(neighbor_system, particle_system),
+                               smoothing_length_particle, sound_speed)
+    nu_b = kinematic_viscosity(neighbor_system,
+                               viscosity_model(particle_system, neighbor_system),
+                               smoothing_length_neighbor, sound_speed)
+
+    v_a = viscous_velocity(v_particle_system, particle_system, particle)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_diff = v_a - v_b
+
+    # ------------------------------------------------------------------------------
+    # SGS part: Compute the subgrid-scale eddy viscosity.
+    # ------------------------------------------------------------------------------
+    # In classical LES the Smagorinsky model defines:
+    #   ν_SGS = (C_S Δ)^2 |S|,
+    # where |S| is the norm of the strain-rate tensor Sᵢⱼ = ½(∂ᵢvⱼ+∂ⱼvᵢ).
+    #
+    # In SPH, one could compute ∂ᵢvⱼ via kernel gradients, but this is costly.
+    # A common low-order surrogate is to approximate the strain‐rate magnitude by a
+    # finite-difference along each particle pair:
+    #
+    #   |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ),
+    #
+    # where δ regularizes the denominator to avoid singularities when particles are very close.
+    #
+    # This yields:
+    #   S_mag = norm(v_diff) / (distance + ε)
+    #
+    # and then the Smagorinsky eddy viscosity:
+    #   ν_SGS = (C_S * h̄)^2 * S_mag.
+    S_mag = norm(v_diff) / (distance + epsilon)
+    nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag
+
+    # Effective viscosities include the SGS term.
+    nu_a_eff = nu_a + nu_SGS
+    nu_b_eff = nu_b + nu_SGS
+
+    # For the Morris model, dynamic viscosities are:
+    mu_a = nu_a_eff * rho_a
+    mu_b = nu_b_eff * rho_b
+
+    force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) /
+                   (distance^2 + epsilon * smoothing_length_average^2) * v_diff
+    return m_b * force_Morris
+end
+
+function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length,
+                             sound_speed)
+    return viscosity.nu
 end