Add simple SGS to Adami and Morris Viscosity

NEWS.md

@@ -8,6 +8,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 ### Features
 
+- 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
   when combined with thread pinning (#722).

docs/src/refs.bib

@@ -715,7 +715,7 @@ @article{MeloEspinosa2014
   journal = {Energy Conversion and Management},
   volume = {84},
   pages = {50--60},
-  doi = {https://doi.org/10.1016/j.enconman.2014.08.004},
+  doi = {10.1016/j.enconman.2014.08.004},
   year = {2014}
 }
 
@@ -736,6 +736,28 @@ @book{Poling2001
   address = {New York}
 }
 
+@article{Smagorinsky1963,
+  author    = {Smagorinsky, Joseph},
+  title     = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment},
+  journal   = {Monthly Weather Review},
+  volume    = {91},
+  number    = {3},
+  pages     = {99--164},
+  year      = {1963},
+  doi       = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
+}
+
+@inproceedings{Lilly1967,
+  author    = {Lilly, Douglas K.},
+  title     = {The Representation of Small‑Scale Turbulence in Numerical Simulation Experiments},
+  booktitle = {Proceedings of the {IBM} Scientific Computing Symposium on Environmental Sciences},
+  editor    = {Goldstine, Herman H.},
+  address   = {Yorktown Heights, New York},
+  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},

examples/fluid/dam_break_2d.jl

@@ -47,18 +47,29 @@ 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)
-# viscosity = ViscosityAdami(nu=nu)
+alpha = 0.02
+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
+
+# Alternatively a kinematic viscosity for water can be set
+# nu = 1.0e-6
+
+# This allows the use of a physical viscosity model like:
+# viscosity_fluid = ViscosityAdami(nu=nu)
+# or with additional dissipation through a Smagorinsky model
+# 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.
 # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
 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,
@@ -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)
 

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)
 

examples/fluid/taylor_green_vortex_2d.jl

@@ -3,7 +3,7 @@
 # P. Ramachandran, K. Puri
 # "Entropically damped artificial compressibility for SPH".
 # In: Computers and Fluids, Volume 179 (2019), pages 579-594.
-# https://doi.org/10.1016/j.compfluid.2018.11.023
+# https://doi.org/10.1016/j.compfluid.2018.11.023
 
 using TrixiParticles
 using OrdinaryDiffEq

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 tensile_instability_control

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