Add simple SGS to Adami and Morris Viscosity

NEWS.md

@@ -8,7 +8,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 ### Features
 
-- New Viscosity model 'ViscosityMorrisSGS' and 'ViscosityAdamiSGS' were added which use simplified Smagorinsky type SGS (#753).
+- New viscosity models `ViscosityMorrisSGS` and `ViscosityAdamiSGS` were added, which use a simplified Smagorinsky-type SGS (#753).  
 
 - With all CPU backends, a new array type is used for the integration array, which defines
   broadcasting to be multithreaded, leading to speedups of up to 5x with large thread counts

docs/src/refs.bib

@@ -756,4 +756,18 @@ @inproceedings{Lilly1967
   pages     = {195--210},
   year      = {1967},
   doi       = {10.5065/D62R3PMM}
-}
+}
+
+@article{Zhu2021,
+  author    = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu},
+  journal   = {Journal of Hydrodynamics},
+  title     = {A {CAD}-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction},
+  year      = {2021},
+  issn      = {1878-0342},
+  month     = apr,
+  number    = {2},
+  pages     = {195--206},
+  volume    = {33},
+  doi       = {10.1007/s42241-021-0031-y},
+  publisher = {Springer Science and Business Media LLC}
+}

examples/fluid/dam_break_2d.jl

@@ -48,7 +48,7 @@ smoothing_kernel = WendlandC2Kernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
 alpha = 0.02
-viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
+viscosity_fluid = 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
@@ -57,10 +57,10 @@ viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
 # nu = 1.0e-6
 
 # This allows the use of a physical viscosity model like:
-# viscosity = ViscosityAdami(nu=nu)
+# viscosity_fluid = 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).
+# viscosity_fluid = ViscosityAdamiSGS(nu=nu)
+# For more details see the documentation "Viscosity model overview".
 
 # Alternatively the density diffusion model by Molteni & Colagrossi can be used,
 # which will run faster.
@@ -69,7 +69,7 @@ density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
-                                           smoothing_length, viscosity=viscosity,
+                                           smoothing_length, viscosity=viscosity_fluid,
                                            density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity), correction=nothing,
                                            surface_tension=nothing,
@@ -80,8 +80,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
 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)
+# 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,

examples/fluid/dam_break_2phase_2d.jl

@@ -37,7 +37,7 @@ water_viscosity = ViscosityMorris(nu=nu_sim_water)
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
-              viscosity=water_viscosity, smoothing_length=smoothing_length,
+              viscosity_fluid=water_viscosity, smoothing_length=smoothing_length,
               gravity=gravity, tspan=tspan, density_diffusion=nothing,
               sound_speed=sound_speed, exponent=7,
               tank_size=(floor(5.366 * H / fluid_particle_spacing) * fluid_particle_spacing,

examples/fluid/dam_break_oil_film_2d.jl

@@ -30,7 +30,8 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 # TODO: broken if both systems use surface tension
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing, tspan=tspan,
-              viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
+              viscosity_fluid=ViscosityMorris(nu=nu_sim_water),
+              smoothing_length=smoothing_length,
               gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed,
               prefix="", reference_particle_spacing=fluid_particle_spacing)
 

examples/fluid/hydrostatic_water_column_2d.jl

@@ -32,15 +32,15 @@ smoothing_length = 1.2 * fluid_particle_spacing
 smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 alpha = 0.02
-viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
+viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
 
 fluid_density_calculator = ContinuityDensity()
 
 # This is to set acceleration with `trixi_include`
 system_acceleration = (0.0, -gravity)
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
-                                           smoothing_length, viscosity=viscosity,
+                                           smoothing_length, viscosity=viscosity_fluid,
                                            acceleration=system_acceleration,
                                            source_terms=nothing)
 

src/schemes/fluid/viscosity.jl

@@ -245,32 +245,32 @@ by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Sm
 The effective kinematic viscosity is defined as
 
 ```math
-\nu_{\mathrm{eff}} = \nu_{\\mathrm{std}} + \nu_{\\mathrm{SGS}},
+\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}},
 ```
 
 with
 
 ```math
-\nu_{\mathrm{SGS}} = (C_S * h)^2 * |S|,
+\nu_{\text{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},
+|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon},
 ```
 
 where:
-- $C_S$ is the Smagorinsky constant (typically 0.1 to 0.2),
-- $h$ is the local smoothing length, and
+- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2),
+- ``h`` is the local smoothing length.
 
 The effective dynamic viscosities are then computed as
 ```math
-\eta_{a,\mathrm{eff}} = \rho_a\, \nu_{\mathrm{eff}},
+\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{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}}}.
+\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}.
 ```
 
 This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation.
@@ -281,12 +281,11 @@ This model is appropriate for turbulent flows where unresolved scales contribute
 - `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]
+    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,
@@ -322,9 +321,9 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps
     #
     # 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:
+    # finite difference along each particle pair:
     #
-    #   |S| ≈ ‖vₐ–v_b‖ / (‖rₐ–r_b‖ + δ),
+    #   |S| ≈ ‖v_ab‖ / (‖r_ab‖ + δ),
     #
     # where δ regularizes the denominator to avoid singularities when particles are very close.
     #
@@ -354,24 +353,36 @@ end
     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 [Smagorinsky (1963)](@cite Smagorinsky1963) for modeling turbulent flows.
+by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure.
+The effective kinematic viscosity is defined as
+
+```math
+\nu_{\text{eff}} = \nu_{\text{std}} + \nu_{\text{SGS}},
+```
+
+with
 
-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.
+\nu_{\text{SGS}} = (C_S h)^2 |S|,
+```
+
+and an approximation for the strain rate magnitude given by
 
-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}})
+|S| \approx \frac{\|v_ab\|}{\|r_ab\| + \epsilon},
 ```
-The standard kinematic viscosity $\nu_{\mathrm{std}}$ is provided by the `nu` parameter. The SGS kinematic viscosity is calculated using the Smagorinsky model:
+
+where:
+- ``C_S`` is the Smagorinsky constant (typically 0.1 to 0.2),
+- ``h`` is the local smoothing length.
+
+The effective dynamic viscosities are then computed as
 ```math
-\nu_{i,\mathrm{SGS}} = (C_S h_{ab})^2 |\bar{S}_{ab}|
+\eta_{a,\text{eff}} = \rho_a\, \nu_{\text{eff}},
 ```
-where $C_S$ is the Smagorinsky constant. This implementation uses a simplified pairwise approximation for the strain rate tensor magnitude $|\bar{S}_{ab}|$:
+and averaged as
 ```math
-|\bar{S}_{ab}| \approx \frac{\|v_a - v_b\|}{\|r_a - r_b\| + \epsilon}
+\bar{\eta}_{ab} = \frac{2 \eta_{a,\text{eff}} \eta_{b,\text{eff}}}{\eta_{a,\text{eff}}+\eta_{b,\text{eff}}}.
 ```
 
 This model is appropriate for turbulent flows where unresolved scales contribute additional dissipation.
@@ -413,26 +424,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e
     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 [Lilly (1967)](@cite Lilly1967) 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.
+    # See comments above for `ViscosityAdamiSGS`.
     S_mag = norm(v_diff) / (distance + epsilon)
     nu_SGS = (viscosity.C_S * smoothing_length_average)^2 * S_mag
 

test/examples/examples_fluid.jl

@@ -25,18 +25,18 @@
                                               clip_negative_pressure=true),
             "WCSPH with ViscosityAdami" => (
                                             # from 0.02*10.0*1.2*0.05/8
-                                            viscosity=ViscosityAdami(nu=0.0015),),
+                                            viscosity_fluid=ViscosityAdami(nu=0.0015),),
             "WCSPH with ViscosityMorris" => (
                                              # from 0.02*10.0*1.2*0.05/8
-                                             viscosity=ViscosityMorris(nu=0.0015),),
+                                             viscosity_fluid=ViscosityMorris(nu=0.0015),),
             "WCSPH with ViscosityAdami and SummationDensity" => (
                                                                  # from 0.02*10.0*1.2*0.05/8
-                                                                 viscosity=ViscosityAdami(nu=0.0015),
+                                                                 viscosity_fluid=ViscosityAdami(nu=0.0015),
                                                                  fluid_density_calculator=SummationDensity(),
                                                                  clip_negative_pressure=true),
             "WCSPH with ViscosityMorris and SummationDensity" => (
                                                                   # from 0.02*10.0*1.2*0.05/8
-                                                                  viscosity=ViscosityMorris(nu=0.0015),
+                                                                  viscosity_fluid=ViscosityMorris(nu=0.0015),
                                                                   fluid_density_calculator=SummationDensity(),
                                                                   clip_negative_pressure=true),
             "WCSPH with smoothing_length=1.3" => (smoothing_length=1.3,),
@@ -55,15 +55,15 @@
                                                                                          smoothing_kernel,
                                                                                          smoothing_length,
                                                                                          sound_speed,
-                                                                                         viscosity=viscosity,
+                                                                                         viscosity=viscosity_fluid,
                                                                                          density_calculator=ContinuityDensity(),
                                                                                          acceleration=(0.0,
                                                                                                        -gravity))),
             "EDAC with SummationDensity" => (fluid_system=EntropicallyDampedSPHSystem(tank.fluid,
                                                                                       smoothing_kernel,
                                                                                       smoothing_length,
                                                                                       sound_speed,
-                                                                                      viscosity=viscosity,
+                                                                                      viscosity=viscosity_fluid,
                                                                                       density_calculator=SummationDensity(),
                                                                                       acceleration=(0.0,
                                                                                                     -gravity)),),

test/schemes/fluid/viscosity.jl

@@ -1,4 +1,3 @@
-include("../../test_util.jl")
 @testset verbose=true "Viscosity" begin
     particle_spacing = 0.2
     smoothing_length = 1.2 * particle_spacing