Fix viscosity calculation when using two different fluids (#569)

* extract changes from #564

* fix

* fix nus

* fix

* add example

* fix comment

* format

* fix tests

* fix tlsph

* fix error

* format

* remove unecessary include

* implement review suggestions

* Fix examples.jl after Case move

* Forgot a part

* fix example

* review implementation

* format & example restruct

* fix test

* fix tests

* fix

* for consitency lets change the nu assignment

* fix tests

* implement suggestion

* change comment

* change comment again

examples/fluid/dam_break_oil_film_2d.jl

@@ -0,0 +1,61 @@
+# 2D dam break simulation with an oil film on top
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# Size parameters
+H = 0.6
+W = 2 * H
+
+gravity = 9.81
+tspan = (0.0, 5.7 / sqrt(gravity))
+
+# Resolution
+fluid_particle_spacing = H / 60
+
+# Numerical settings
+smoothing_length = 3.5 * fluid_particle_spacing
+sound_speed = 20 * sqrt(gravity * H)
+
+nu = 0.02 * smoothing_length * sound_speed / 8
+oil_viscosity = ViscosityMorris(nu=10 * nu)
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
+              sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
+              viscosity=ViscosityMorris(nu=nu), smoothing_length=smoothing_length,
+              gravity=gravity,
+              density_diffusion=DensityDiffusionMolteniColagrossi(delta=0.1),
+              sound_speed=sound_speed)
+
+# ==========================================================================================
+# ==== Setup oil layer
+
+oil_size = (W, 0.1 * H)
+oil_density = 700.0
+
+oil = RectangularShape(fluid_particle_spacing,
+                       round.(Int, oil_size ./ fluid_particle_spacing),
+                       zeros(length(oil_size)), density=oil_density)
+
+# move on top of the water
+for i in axes(oil.coordinates, 2)
+    oil.coordinates[:, i] .+= [0.0, H]
+end
+
+oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
+                                         state_equation, smoothing_kernel,
+                                         smoothing_length, viscosity=oil_viscosity,
+                                         density_diffusion=density_diffusion,
+                                         acceleration=(0.0, -gravity),
+                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02),
+                                         correction=AkinciFreeSurfaceCorrection(oil_density))
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(fluid_system, oil_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch{2}(update_strategy=nothing))
+ode = semidiscretize(semi, tspan, data_type=nothing)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # This is overwritten by the stepsize callback
+            save_everystep=false, callback=callbacks);

examples/fluid/falling_water_spheres_2d.jl

@@ -90,6 +90,6 @@ callbacks = CallbackSet(info_callback, saving_callback)
 
 # Use a Runge-Kutta method with automatic (error based) time step size control.
 sol = solve(ode, RDPK3SpFSAL35(),
-            abstol=1e-6, # Default abstol is 1e-6
+            abstol=1e-7, # Default abstol is 1e-6
             reltol=1e-4, # Default reltol is 1e-3
             save_everystep=false, callback=callbacks);

src/schemes/boundary/boundary.jl

@@ -1,4 +1,6 @@
 include("dummy_particles/dummy_particles.jl")
 include("system.jl")
+include("open_boundary/boundary_zones.jl")
+include("open_boundary/system.jl")
 # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`
 # and the `TotalLagrangianSPHSystem` and are therefore included later.

src/schemes/boundary/open_boundary/system.jl

@@ -481,3 +481,10 @@ end
 function wrap_reference_function(constant_vector_, ::Val{NDIMS}) where {NDIMS}
     return constant_vector(coords, t) = SVector{NDIMS}(constant_vector_)
 end
+
+# To account for boundary effects in the viscosity term of the RHS, use the viscosity model
+# of the neighboring particle systems.
+@inline viscosity_model(system::OpenBoundarySPHSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity
+@inline viscosity_model(system::OpenBoundarySPHSystem, neighbor_system::BoundarySystem) = neighbor_system.boundary_model.viscosity
+# When the neighbor is an open boundary system, just use the viscosity of the fluid `system` instead
+@inline viscosity_model(system, neighbor_system::OpenBoundarySPHSystem) = system.viscosity

src/schemes/boundary/system.jl

@@ -383,9 +383,9 @@ function restart_with!(system::BoundarySPHSystem{<:BoundaryModelDummyParticles{C
     return system
 end
 
-function viscosity_model(system::BoundarySPHSystem)
-    return system.boundary_model.viscosity
-end
+# To incorporate the effect at boundaries in the viscosity term of the RHS the neighbor
+# viscosity model has to be used.
+@inline viscosity_model(system::BoundarySPHSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity
 
 function calculate_dt(v_ode, u_ode, cfl_number, system::BoundarySystem)
     return Inf

src/schemes/fluid/fluid.jl

@@ -37,7 +37,10 @@ end
 
 write_v0!(v0, system, density_calculator) = v0
 
-@inline viscosity_model(system::FluidSystem) = system.viscosity
+# To account for boundary effects in the viscosity term of the RHS, use the viscosity model
+# of the neighboring particle systems.
+@inline viscosity_model(system::FluidSystem, neighbor_system::FluidSystem) = neighbor_system.viscosity
+@inline viscosity_model(system::FluidSystem, neighbor_system::BoundarySystem) = neighbor_system.boundary_model.viscosity
 
 function compute_density!(system, u, u_ode, semi, ::ContinuityDensity)
     # No density update with `ContinuityDensity`

src/schemes/fluid/viscosity.jl

@@ -4,20 +4,7 @@ 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_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_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_a, rho_b, grad_kernel)
-    # No viscosity in the open boundary system. Use viscosity of the fluid system.
-    viscosity = viscosity_model(particle_system)
+    viscosity = viscosity_model(particle_system, neighbor_system)
 
     return dv_viscosity(viscosity, particle_system, neighbor_system,
                         v_particle_system, v_neighbor_system,
@@ -154,15 +141,20 @@ end
     v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
     v_diff = v_a - v_b
 
+    nu_a = kinematic_viscosity(particle_system,
+                               viscosity_model(neighbor_system, particle_system))
+    nu_b = kinematic_viscosity(neighbor_system,
+                               viscosity_model(particle_system, neighbor_system))
+
     pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
-                      smoothing_length, grad_kernel)
+                      smoothing_length, grad_kernel, nu_a, nu_b)
 
     return m_b * pi_ab
 end
 
 @inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance,
                                                           rho_mean, rho_a, rho_b, h,
-                                                          grad_kernel)
+                                                          grad_kernel, nu_a, nu_b)
     (; alpha, beta, epsilon) = viscosity
 
     # v_ab ⋅ r_ab
@@ -181,12 +173,12 @@ end
 end
 
 @inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean,
-                                              rho_a, rho_b, h, grad_kernel)
-    (; epsilon, nu) = viscosity
+                                              rho_a, rho_b, h, grad_kernel, nu_a,
+                                              nu_b)
+    epsilon = viscosity.epsilon
 
-    # 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
+    mu_a = nu_a * rho_a
+    mu_b = nu_b * rho_b
 
     return (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) /
            (distance^2 + epsilon * h^2) * v_diff
@@ -247,16 +239,20 @@ end
                                              particle, neighbor, pos_diff,
                                              distance, sound_speed, m_a, m_b,
                                              rho_a, rho_b, grad_kernel)
-    (; epsilon, nu) = viscosity
     (; smoothing_length) = particle_system
 
+    epsilon = viscosity.epsilon
+    nu_a = kinematic_viscosity(particle_system,
+                               viscosity_model(neighbor_system, particle_system))
+    nu_b = kinematic_viscosity(neighbor_system,
+                               viscosity_model(particle_system, neighbor_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
 
-    # 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
+    eta_a = nu_a * rho_a
+    eta_b = nu_b * rho_b
 
     eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b)
 

src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl

@@ -51,7 +51,8 @@ struct DensityDiffusionMolteniColagrossi{ELTYPE} <: DensityDiffusion
 end
 
 @inline function density_diffusion_psi(::DensityDiffusionMolteniColagrossi, rho_a, rho_b,
-                                       pos_diff, distance, system, particle, neighbor)
+                                       pos_diff, distance, system, neighbor_system,
+                                       particle, neighbor)
     return 2 * (rho_a - rho_b) * pos_diff / distance^2
 end
 
@@ -86,7 +87,8 @@ struct DensityDiffusionFerrari <: DensityDiffusion
 end
 
 @inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b,
-                                       pos_diff, distance, system, particle, neighbor)
+                                       pos_diff, distance, system, neighbor_system,
+                                       particle, neighbor)
     (; smoothing_length) = system
 
     return ((rho_a - rho_b) / 2smoothing_length) * pos_diff / distance
@@ -170,12 +172,13 @@ function Base.show(io::IO, density_diffusion::DensityDiffusionAntuono)
 end
 
 @inline function density_diffusion_psi(density_diffusion::DensityDiffusionAntuono,
-                                       rho_a, rho_b,
-                                       pos_diff, distance, system, particle, neighbor)
+                                       rho_a, rho_b, pos_diff, distance, system,
+                                       neighbor_system, particle, neighbor)
     (; normalized_density_gradient) = density_diffusion
 
     normalized_gradient_a = extract_svector(normalized_density_gradient, system, particle)
-    normalized_gradient_b = extract_svector(normalized_density_gradient, system, neighbor)
+    normalized_gradient_b = extract_svector(normalized_density_gradient, neighbor_system,
+                                            neighbor)
 
     # Fist term by Molteni & Colagrossi
     result = 2 * (rho_a - rho_b)
@@ -242,7 +245,7 @@ end
     volume_b = m_b / rho_b
 
     psi = density_diffusion_psi(density_diffusion, rho_a, rho_b, pos_diff, distance,
-                                particle_system, particle, neighbor)
+                                particle_system, neighbor_system, particle, neighbor)
     density_diffusion_term = dot(psi, grad_kernel) * volume_b
 
     dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term

src/schemes/schemes.jl

@@ -1,8 +1,5 @@
 # Include all schemes without rhs first. The rhs depends on the systems to define
 # interactions between the different system types.
-# Viscosity requires the open boundary system.
-include("boundary/open_boundary/boundary_zones.jl")
-include("boundary/open_boundary/system.jl")
 include("fluid/fluid.jl")
 include("boundary/boundary.jl")
 include("solid/total_lagrangian_sph/total_lagrangian_sph.jl")

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -379,10 +379,6 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u)
     restart_with!(system, system.boundary_model, v, u)
 end
 
-function viscosity_model(system::TotalLagrangianSPHSystem)
-    return system.boundary_model.viscosity
-end
-
 # An explanation of these equation can be found in
 # J. Lubliner, 2008. Plasticity theory.
 # See here below Equation 5.3.21 for the equation for the equivalent stress.
@@ -434,3 +430,8 @@ function cauchy_stress(system::TotalLagrangianSPHSystem)
 
     return cauchy_stress_tensors
 end
+
+# To account for boundary effects in the viscosity term of the RHS, use the viscosity model
+# of the neighboring particle systems.
+@inline viscosity_model(system::TotalLagrangianSPHSystem, neighbor_system) = neighbor_system.viscosity
+@inline viscosity_model(system::FluidSystem, neighbor_system::TotalLagrangianSPHSystem) = neighbor_system.boundary_model.viscosity

test/examples/examples.jl

@@ -87,6 +87,18 @@
             @test count_rhs_allocations(sol, semi) == 0
         end
 
+        @trixi_testset "fluid/dam_break_oil_film_2d.jl" begin
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "dam_break_oil_film_2d.jl"),
+                                           tspan=(0.0, 0.1)) [
+                r"┌ Info: The desired tank length in y-direction .*\n",
+                r"└ New tank length in y-direction.*\n",
+            ]
+            @test sol.retcode == ReturnCode.Success
+            @test count_rhs_allocations(sol, semi) == 0
+        end
+
         @trixi_testset "fluid/dam_break_2d.jl with KernelAbstractions.jl" begin
             # Emulate the GPU code on the CPU by passing `data_type = Array`
             @test_nowarn_mod trixi_include(@__MODULE__,

test/schemes/fluid/viscosity.jl

@@ -9,42 +9,54 @@
 
     fluid = rectangular_patch(particle_spacing, (3, 3), seed=1)
 
-    system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
-                                               state_equation, smoothing_kernel,
-                                               smoothing_length)
-
     v_diff = [0.1, -0.75]
     pos_diff = [-0.5 * smoothing_length, 0.75 * smoothing_length]
     distance = norm(pos_diff)
     rho_a = rho_b = rho_mean = 1000.0
 
-    grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
-                                                       distance)
-
     # We only test here that the values don't change
     @testset verbose=true "`ArtificialViscosityMonaghan`" begin
         viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
 
+        system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length, viscosity=viscosity)
+
+        grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
+                                                           distance)
+
         dv = viscosity(sound_speed, v_diff, pos_diff, distance,
                        rho_mean, rho_a, rho_b, smoothing_length,
-                       grad_kernel)
+                       grad_kernel, 0.0, 0.0)
 
         @test isapprox(dv[1], -0.007073849138494646, atol=6e-15)
         @test isapprox(dv[2], 0.01061077370774197, atol=6e-15)
     end
     @testset verbose=true "`ViscosityMorris`" begin
-        viscosity = ViscosityMorris(nu=7e-3)
+        nu = 7e-3
+        viscosity = ViscosityMorris(nu=nu)
+        system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length, viscosity=viscosity)
+
+        grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
+                                                           distance)
 
         dv = viscosity(sound_speed, v_diff, pos_diff, distance,
                        rho_mean, rho_a, rho_b, smoothing_length,
-                       grad_kernel)
+                       grad_kernel, nu, nu)
 
         @test isapprox(dv[1], -0.00018421886647594437, atol=6e-15)
         @test isapprox(dv[2], 0.0013816414985695826, atol=6e-15)
     end
     @testset verbose=true "`ViscosityAdami`" begin
         viscosity = ViscosityAdami(nu=7e-3)
+        system_wcsph = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length, viscosity=viscosity)
 
+        grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
+                                                           distance)
         v = fluid.velocity
 
         m_a = 0.01