Add Viscosity Formulation as per Morris et al.

src/TrixiParticles.jl

@@ -52,7 +52,7 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel
 export StateEquationCole
-export ArtificialViscosityMonaghan, ViscosityAdami
+export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMoris
 export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
        DensityDiffusionAntuono
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,

src/schemes/fluid/viscosity.jl

@@ -29,17 +29,17 @@ function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system,
 end
 
 @doc raw"""
-    ArtificialViscosityMonaghan(; alpha, beta, epsilon=0.01)
+    ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01)
 
 # Keywords
 - `alpha`: A value of `0.02` is usually used for most simulations. For a relation with the
            kinematic viscosity, see description below.
-- `beta`: A value of `0.0` works well for simulations with shocks of moderate strength.
-          In simulations where the Mach number can be very high, eg. astrophysical calculation,
-          good results can be obtained by choosing a value of `beta=2` and `alpha=1`.
+- `beta=0.0`: A value of `0.0` works well for most fluid simulations and simulations
+          with shocks of moderate strength. In simulations where the Mach number can be
+          very high, eg. astrophysical calculation, good results can be obtained by
+          choosing a value of `beta=2` and `alpha=1`.
 - `epsilon=0.01`: Parameter to prevent singularities.
 
-
 Artificial viscosity by Monaghan (Monaghan 1992, Monaghan 1989), given by
 ```math
 \Pi_{ab} =
@@ -80,7 +80,7 @@ struct ArtificialViscosityMonaghan{ELTYPE}
     beta    :: ELTYPE
     epsilon :: ELTYPE
 
-    function ArtificialViscosityMonaghan(; alpha, beta, epsilon=0.01)
+    function ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01)
         new{typeof(alpha)}(alpha, beta, epsilon)
     end
 end
@@ -97,7 +97,10 @@ end
     v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
     v_diff = v_a - v_b
 
-    pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, smoothing_length)
+    # v_ab ⋅ r_ab
+    vr = dot(v_diff, pos_diff)
+
+    pi_ab = viscosity(sound_speed, vr, distance, rho_mean, smoothing_length)
 
     if pi_ab < eps()
         return zero(pos_diff)
@@ -106,13 +109,9 @@ end
     return -m_b * pi_ab * smoothing_kernel_grad(particle_system, pos_diff, distance)
 end
 
-@inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance,
-                                                          rho_mean, h)
+@inline function (viscosity::ArtificialViscosityMonaghan)(c, vr, distance, rho_mean, h)
     (; alpha, beta, epsilon) = viscosity
 
-    # v_ab ⋅ r_ab
-    vr = dot(v_diff, pos_diff)
-
     # Monaghan 2005 p. 1741 (doi: 10.1088/0034-4885/68/8/r01):
     # "In the case of shock tube problems, it is usual to turn the viscosity on for
     # approaching particles and turn it off for receding particles. In this way, the
@@ -126,10 +125,9 @@ end
 end
 
 # See, e.g.,
-# M. Antuono, A. Colagrossi, S. Marrone.
-# "Numerical Diffusive Terms in Weakly-Compressible SPH Schemes."
-# In: Computer Physics Communications 183, no. 12 (2012), pages 2570-80.
-# https://doi.org/10.1016/j.cpc.2012.07.006
+# Joseph J. Monaghan. "Smoothed Particle Hydrodynamics".
+# In: Reports on Progress in Physics (2005), pages 1703-1759.
+# [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01)
 function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan)
     (; smoothing_length) = system
     (; alpha) = viscosity
@@ -190,6 +188,7 @@ end
     rho_a = particle_density(v_particle_system, particle_system, particle)
     rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
 
+
     # TODO This is not correct for two different fluids. It should be `nu_a` and `nu_b`.
     eta_a = nu * rho_a
     eta_b = nu * rho_b
@@ -222,4 +221,66 @@ function kinematic_viscosity(system, viscosity::ViscosityAdami)
     return viscosity.nu
 end
 
+@doc raw"""
+    ViscosityMoris(; nu, epsilon=0.01)
+
+Viscosity by Moris (Moris et al. 1997).
+```math
+f_{ab} = \sum_w \frac{m_b(eta_a+eta_b)}{||r_{ab}||^2+(\epsilon h_{ab})^2}  \nabla W_{ab} \cdot r_{ab}\cdot v_{ab},
+```
+where ``\eta_a = \rho_a \nu_a`` with ``\nu`` as the kinematic viscosity.
+
+# Keywords
+- `nu`: Kinematic viscosity
+- `epsilon=0.01`: Parameter to prevent singularities
+
+## References
+- J. Morris et al., "Modeling Low Reynolds Number Incompressible Flows Using SPH",
+  In: Journal of Computational Physics, Volume 136, Issue 1, 1997, Pages 214-226.
+  [doi: doi.org/10.1006/jcph.1997.5776](https://doi.org/10.1006/jcph.1997.5776)
+- G. Fourtakas et al., "Local uniform stencil (LUST) boundary condition for arbitrary
+  3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models",
+  In: Computers & Fluids, 2019.
+  [doi: doi.org/10.1016/j.compfluid.2019.06.009](https://doi.org/10.1016/j.compfluid.2019.06.009)
+"""
+struct ViscosityMoris{ELTYPE}
+    nu::ELTYPE
+    epsilon::ELTYPE
+
+    function ViscosityMoris(; nu, epsilon=0.01)
+        new{typeof(nu)}(nu, epsilon)
+    end
+end
+
+@inline function (viscosity::ViscosityMoris)(particle_system, neighbor_system,
+                                             v_particle_system, v_neighbor_system,
+                                             particle, neighbor, pos_diff,
+                                             distance, sound_speed, m_a, m_b, rho_mean)
+    (; epsilon, nu) = viscosity
+    (; smoothing_length) = particle_system
+
+    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
+
+    rho_a = particle_density(v_particle_system, particle_system, particle)
+    rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
+
+    # TODO This is not correct for two different fluids this should be nu_a and nu_b
+    eta_a = nu * rho_a
+    eta_b = nu * rho_b
+
+    factor = (m_b * (eta_a + eta_b)) /
+             (rho_b * (distance^2 + (epsilon * smoothing_length)^2))
+
+    grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance)
+    visc = factor * dot(pos_diff, grad_kernel) .* v_diff
+
+    return visc
+end
+
+function kinematic_viscosity(system, viscosity::ViscosityMoris)
+    return viscosity.nu
+end
+
 @inline viscous_velocity(v, system, particle) = current_velocity(v, system, particle)

src/visualization/write2vtk.jl

@@ -218,7 +218,7 @@ end
 
 write2vtk!(vtk, viscosity::Nothing) = vtk
 
-function write2vtk!(vtk, viscosity::ViscosityAdami)
+function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMoris})
     vtk["viscosity_nu"] = viscosity.nu
     vtk["viscosity_epsilon"] = viscosity.epsilon
 end
@@ -305,14 +305,17 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, viscosity,
 end
 
 function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles,
-                    viscosity::ViscosityAdami, system; write_meta_data=true)
+                    viscosity::Union{ViscosityAdami, ViscosityMoris}, system;
+                    write_meta_data=true)
     vtk["hydrodynamic_density"] = [particle_density(v, system, particle)
                                    for particle in eachparticle(system)]
     vtk["pressure"] = model.pressure
     vtk["wall_velocity"] = view(model.cache.wall_velocity, 1:ndims(system), :)
 
-    if write_meta_data
+    if write_meta_data && viscosity isa ViscosityAdami
         vtk["viscosity_model"] = "ViscosityAdami"
+    else
+        vtk["viscosity_model"] = "ViscosityMoris"
     end
 
     return vtk