Add Viscosity Formulation as per Morris et al. (#504)

* Add Viscosity Moris * format * fix * implement suggestions * other suggestions * fix equation * format * review comments * fixes * fix * update test * update * fix open boundary * update reference files * update validation criteria * update reference file * new reference files * ups * update ref file * update validation error bounds * update error bound * update error bound * add comment and correct setting of 'v'
trixi-framework · Jun 28, 2024 · ffe3663 · ffe3663
1 parent 708ad9d
commit ffe3663
Show file tree

Hide file tree

Showing 16 changed files with 870 additions and 752 deletions.
diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
@@ -48,6 +48,9 @@ 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)
 # Alternatively the density diffusion model by Molteni & Colagrossi can be used,
 # which will run faster.
 # density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)

diff --git a/examples/fluid/oscillating_drop_2d.jl b/examples/fluid/oscillating_drop_2d.jl
@@ -89,5 +89,5 @@ exact_solution_ode = ODEProblem(exact_solution_rhs, exact_u0, tspan)
 sol_exact = solve(exact_solution_ode, RDPK3SpFSAL49(), save_everystep=false)
 
 # Error in the semi-major axis of the elliptical drop
-error_A = maximum(sol.u[end].x[2]) + 0.5fluid_particle_spacing -
+error_A = maximum(sol.u[end].x[2]) + 0.5 * fluid_particle_spacing -
           maximum(sol_exact.u[end][2:3])
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -58,7 +58,7 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel
 export StateEquationCole
-export ArtificialViscosityMonaghan, ViscosityAdami
+export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris
 export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari,
        DensityDiffusionAntuono
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,

diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -17,7 +17,6 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         rho_a = particle_density(v_particle_system, particle_system, particle)
         rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
-        rho_mean = 0.5 * (rho_a + rho_b)
 
         p_a = particle_pressure(v_particle_system, particle_system, particle)
         p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)
@@ -34,7 +33,7 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_viscosity_ = dv_viscosity(particle_system, neighbor_system,
                                      v_particle_system, v_neighbor_system,
                                      particle, neighbor, pos_diff, distance,
-                                     sound_speed, m_a, m_b, rho_mean)
+                                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 
         for i in 1:ndims(particle_system)
             dv[i, particle] += dv_pressure[i] + dv_viscosity_[i]

diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl
@@ -1,57 +1,49 @@
 
+# Unpack the neighboring systems viscosity to dispatch on the viscosity type
 function dv_viscosity(particle_system, neighbor_system,
                       v_particle_system, v_neighbor_system,
                       particle, neighbor, pos_diff, distance,
-                      sound_speed, m_a, m_b, rho_mean)
+                      sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
     viscosity = viscosity_model(neighbor_system)
 
     return dv_viscosity(viscosity, particle_system, neighbor_system,
                         v_particle_system, v_neighbor_system,
                         particle, neighbor, pos_diff, distance,
-                        sound_speed, m_a, m_b, rho_mean)
+                        sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 end
 
 function dv_viscosity(particle_system, neighbor_system::OpenBoundarySPHSystem,
                       v_particle_system, v_neighbor_system,
                       particle, neighbor, pos_diff, distance,
-                      sound_speed, m_a, m_b, rho_mean)
+                      sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
     # No viscosity in the open boundary system. Use viscosity of the fluid system.
     viscosity = viscosity_model(particle_system)
 
     return dv_viscosity(viscosity, particle_system, neighbor_system,
                         v_particle_system, v_neighbor_system,
                         particle, neighbor, pos_diff, distance,
-                        sound_speed, m_a, m_b, rho_mean)
+                        sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 end
 
 function dv_viscosity(viscosity, particle_system, neighbor_system,
                       v_particle_system, v_neighbor_system,
                       particle, neighbor, pos_diff, distance,
-                      sound_speed, m_a, m_b, rho_mean)
+                      sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
     return viscosity(particle_system, neighbor_system,
                      v_particle_system, v_neighbor_system,
                      particle, neighbor, pos_diff, distance,
-                     sound_speed, m_a, m_b, rho_mean)
+                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 end
 
 function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system,
                       v_particle_system, v_neighbor_system,
                       particle, neighbor, pos_diff, distance,
-                      sound_speed, m_a, m_b, rho_mean)
+                      sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
     return zero(pos_diff)
 end
 
 @doc raw"""
-    ArtificialViscosityMonaghan(; alpha, beta, 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`.
-- `epsilon=0.01`: Parameter to prevent singularities.
-
+    ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01)
 
 Artificial viscosity by Monaghan (Monaghan 1992, Monaghan 1989), given by
 ```math
@@ -77,6 +69,15 @@ To do so, Monaghan (Monaghan 2005) defined an equivalent effective physical kine
 ```
 where ``d`` is the dimension.
 
+# Keywords
+- `alpha`: A value of `0.02` is usually used for most simulations. For a relation with the
+           kinematic viscosity, see description above.
+- `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.0` and `alpha=1.0`.
+- `epsilon=0.01`: Parameter to prevent singularities.
+
 ## References
 - Joseph J. Monaghan. "Smoothed Particle Hydrodynamics".
   In: Annual Review of Astronomy and Astrophysics 30.1 (1992), pages 543-574.
@@ -93,34 +94,75 @@ 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
 
-@inline function (viscosity::ArtificialViscosityMonaghan)(particle_system, neighbor_system,
-                                                          v_particle_system,
-                                                          v_neighbor_system,
-                                                          particle, neighbor, pos_diff,
-                                                          distance, sound_speed, m_a, m_b,
-                                                          rho_mean)
+@doc raw"""
+    ViscosityMorris(; nu, epsilon=0.01)
+
+Viscosity by Morris et al. (1997).
+
+To the force ``f_{ab}`` between two particles ``a`` and ``b`` due to pressure gradients,
+an additional force term ``\tilde{f}_{ab}`` is added with
+```math
+\tilde{f}_{ab} = m_a m_b \frac{(\mu_a + \mu_b) r_{ab} \cdot \nabla W_{ab}}{\rho_a \rho_b (\Vert r_{ab} \Vert^2 + \epsilon h^2)} v_{ab},
+```
+where ``\mu_a = \rho_a \nu`` and ``\mu_b = \rho_b \nu`` denote the dynamic viscosity
+of particle ``a`` and ``b`` respectively, and ``\nu`` is the kinematic viscosity.
+
+# Keywords
+- `nu`: Kinematic viscosity
+- `epsilon=0.01`: Parameter to prevent singularities
+
+## References
+- Joseph P. Morris, Patrick J. Fox, Yi Zhu.
+  "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)
+- Georgios Fourtakas, Jose M. Dominguez, Renato Vacondio, Benedict D. Rogers.
+  "Local uniform stencil (LUST) boundary condition for arbitrary
+  3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models".
+  In: Computers & Fluids, Volume 190 (2019), pages 346--361.
+  [doi: 10.1016/j.compfluid.2019.06.009](https://doi.org/10.1016/j.compfluid.2019.06.009)
+"""
+struct ViscosityMorris{ELTYPE}
+    nu::ELTYPE
+    epsilon::ELTYPE
+
+    function ViscosityMorris(; nu, epsilon=0.01)
+        new{typeof(nu)}(nu, epsilon)
+    end
+end
+
+function kinematic_viscosity(system, viscosity::ViscosityMorris)
+    return viscosity.nu
+end
+
+@inline function (viscosity::Union{ArtificialViscosityMonaghan,
+                                   ViscosityMorris})(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)
     (; smoothing_length) = particle_system
 
+    rho_mean = 0.5 * (rho_a + rho_b)
+
     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
 
-    pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, smoothing_length)
+    pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
+                      smoothing_length, grad_kernel)
 
-    if pi_ab < eps()
-        return zero(pos_diff)
-    end
-
-    return -m_b * pi_ab * smoothing_kernel_grad(particle_system, pos_diff, distance)
+    return m_b * pi_ab
 end
 
 @inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance,
-                                                          rho_mean, h)
+                                                          rho_mean, rho_a, rho_b, h,
+                                                          grad_kernel)
     (; alpha, beta, epsilon) = viscosity
 
     # v_ab ⋅ r_ab
@@ -132,17 +174,28 @@ end
     # viscosity is used for shocks and not rarefactions."
     if vr < 0
         mu = h * vr / (distance^2 + epsilon * h^2)
-        return -(alpha * c * mu + beta * mu^2) / rho_mean
+        return (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel
     end
 
-    return 0.0
+    return zero(v_diff)
+end
+
+@inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean,
+                                              rho_a, rho_b, h, grad_kernel)
+    (; epsilon, nu) = viscosity
+
+    # TODO This is not correct for two different fluids. It should be `nu_a` and `nu_b`.
+    mu_a = nu * rho_a
+    mu_b = nu * rho_b
+
+    return (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) /
+           (distance^2 + epsilon * h^2) * v_diff
 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
@@ -192,17 +245,15 @@ end
 @inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system,
                                              v_particle_system, v_neighbor_system,
                                              particle, neighbor, pos_diff,
-                                             distance, sound_speed, m_a, m_b, rho_mean)
+                                             distance, sound_speed, m_a, m_b,
+                                             rho_a, rho_b, grad_kernel)
     (; 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. It should be `nu_a` and `nu_b`.
     eta_a = nu * rho_a
     eta_b = nu * rho_b
@@ -215,8 +266,6 @@ end
     volume_a = m_a / rho_a
     volume_b = m_b / rho_b
 
-    grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance)
-
     # 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.

diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -56,7 +56,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                         dv_viscosity(particle_system, neighbor_system,
                                      v_particle_system, v_neighbor_system,
                                      particle, neighbor, pos_diff, distance,
-                                     sound_speed, m_a, m_b, rho_mean)
+                                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 
         dv_surface_tension = surface_tension_correction *
                              surface_tension_force(surface_tension_a, surface_tension_b,

diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl
@@ -109,7 +109,7 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_viscosity_ = dv_viscosity(neighbor_system, particle_system,
                                      v_neighbor_system, v_particle_system,
                                      neighbor, particle, pos_diff, distance,
-                                     sound_speed, m_b, m_a, rho_mean)
+                                     sound_speed, m_b, m_a, rho_a, rho_b, grad_kernel)
 
         dv_particle = dv_boundary + dv_viscosity_
 

diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl
@@ -255,7 +255,7 @@ end
 
 write2vtk!(vtk, viscosity::Nothing) = vtk
 
-function write2vtk!(vtk, viscosity::ViscosityAdami)
+function write2vtk!(vtk, viscosity::Union{ViscosityAdami, ViscosityMorris})
     vtk["viscosity_nu"] = viscosity.nu
     vtk["viscosity_epsilon"] = viscosity.epsilon
 end
@@ -373,14 +373,15 @@ 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::ViscosityAdami, 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
-        vtk["viscosity_model"] = "ViscosityAdami"
+        vtk["viscosity_model"] = type2string(viscosity)
     end
 
     return vtk

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
@@ -8,7 +8,7 @@
                                                     "oscillating_drop_2d.jl"))
             @test sol.retcode == ReturnCode.Success
             # This error varies between serial and multithreaded runs
-            @test isapprox(error_A, 0.0001717690010767381, atol=5e-7)
+            @test isapprox(error_A, 0.0, atol=1.73e-4)
             @test count_rhs_allocations(sol, semi) == 0
         end
 

diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl
@@ -2,3 +2,4 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl")
 include("rhs.jl")
 include("pressure_acceleration.jl")
 include("surface_tension.jl")
+include("viscosity.jl")