From cf28d2e6a0c7982bfa9c4566e8fff82dda37f25c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 1 Aug 2024 07:55:31 +0200 Subject: [PATCH] inter --- .../fluid/falling_water_spheres_morris_2d.jl | 94 +++++++++++++ src/TrixiParticles.jl | 2 +- .../fluid/entropically_damped_sph/rhs.jl | 3 +- .../fluid/entropically_damped_sph/system.jl | 19 +-- src/schemes/fluid/surface_normal_sph.jl | 131 +++++++++++++++--- src/schemes/fluid/surface_tension.jl | 52 +++++-- .../fluid/weakly_compressible_sph/rhs.jl | 3 +- .../fluid/weakly_compressible_sph/system.jl | 21 +-- src/visualization/write2vtk.jl | 4 + 9 files changed, 275 insertions(+), 54 deletions(-) create mode 100644 examples/fluid/falling_water_spheres_morris_2d.jl diff --git a/examples/fluid/falling_water_spheres_morris_2d.jl b/examples/fluid/falling_water_spheres_morris_2d.jl new file mode 100644 index 000000000..85a08145f --- /dev/null +++ b/examples/fluid/falling_water_spheres_morris_2d.jl @@ -0,0 +1,94 @@ +# In this example two circles of water drop to the floor demonstrating the difference +# between the behavior with and without surface tension modelling. +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +boundary_layers = 3 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +gravity = 9.81 +tspan = (0.0, 0.3) + +# Boundary geometry and initial fluid particle positions +initial_fluid_size = (0.0, 0.0) +tank_size = (2.0, 0.5) + +fluid_density = 1000.0 +sound_speed = 100 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation=state_equation) + +sphere_radius = 0.05 + +sphere1_center = (0.5, 0.2) +sphere2_center = (1.5, 0.2) +sphere1 = SphereShape(fluid_particle_spacing, sphere_radius, sphere1_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) +sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, + fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, -3.0)) + +# ========================================================================================== +# ==== Fluid +fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps() +fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}() + +fluid_density_calculator = ContinuityDensity() + +nu = 0.005 +alpha = 8 * nu / (fluid_smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0) +density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1) + +sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel, + fluid_smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=ContinuityDensity(), + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.05)) + +sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator, + state_equation, fluid_smoothing_kernel, + fluid_smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +wall_viscosity = nu +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + fluid_smoothing_kernel, fluid_smoothing_length, + viscosity=ViscosityAdami(nu=wall_viscosity)) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, + adhesion_coefficient=1.0) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=50) +saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", + prefix="", write_meta_data=true) + +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-7, # Default abstol is 1e-6 + reltol=1e-4, # Default reltol is 1e-3 + save_everystep=false, callback=callbacks); diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 00b583ac4..14e5d54cd 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -79,7 +79,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk -export SurfaceTensionAkinci, CohesionForceAkinci +export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris export ColorfieldSurfaceNormal end # module diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 42465cc70..361a43e8f 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -40,7 +40,8 @@ function interact!(dv, v_particle_system, u_particle_system, dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, - particle, neighbor, pos_diff, distance) + particle, neighbor, pos_diff, distance, + rho_a, rho_b) dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 20fd48f77..ece424593 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -103,7 +103,10 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, cache = create_cache_density(initial_condition, density_calculator) cache = (; - create_cache_edac(surface_tension, ELTYPE, NDIMS, n_particles)..., + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., cache...) new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), @@ -119,16 +122,6 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, end end -function create_cache_edac(surface_tension_model, ELTYPE, NDIMS, nparticles) - return (;) -end - -function create_cache_edac(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) - surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) - neighbor_count = Array{ELTYPE, 1}(undef, nparticles) - return (; surface_normal, neighbor_count) -end - function Base.show(io::IO, system::EntropicallyDampedSPHSystem) @nospecialize system # reduce precompilation time @@ -205,7 +198,9 @@ end function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) compute_density!(system, u, u_ode, semi, system.density_calculator) - compute_surface_normal!(system, system.surface_tension, v, u, v_ode, u_ode, semi, t) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, + t) + compute_curvature!(system, system.surface_tension, v, u, v_ode, u_ode, semi, t) end function write_v0!(v0, system::EntropicallyDampedSPHSystem, density_calculator) diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index acdb3379f..6fa80019a 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -7,11 +7,26 @@ function ColorfieldSurfaceNormal(; smoothing_kernel, smoothing_length) return ColorfieldSurfaceNormal(smoothing_kernel, smoothing_length) end +function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) + surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_normal, neighbor_count) +end + +function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surfn) + # Normal not needed + return system +end + # Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: This is the simplest form of normal approximation commonly used in SPH and comes -# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries -function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v, - v_neighbor_system, u_neighbor_system, semi, surfn) +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn) (; cache) = system (; smoothing_kernel, smoothing_length) = surfn @@ -37,9 +52,8 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v, grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) for i in 1:ndims(system) # The `smoothing_length` here is used for scaling - # TODO move this to the surface tension model since this is not a general thing cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] * smoothing_length + grad_kernel[i] end cache.neighbor_count[particle] += 1 @@ -49,10 +63,9 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, u_system, v, end # Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: This is the simplest form of normal approximation commonly used in SPH and comes -# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries -function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system, v, - v_neighbor_system, u_neighbor_system, semi, surfn) +# and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_normal!(system::FluidSystem, neighbor_system::BoundarySystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn) (; cache) = system (; colorfield, colorfield_bnd) = neighbor_system.boundary_model.cache (; smoothing_kernel, smoothing_length) = surfn @@ -101,9 +114,8 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system, smoothing_length) for i in 1:ndims(system) # The `smoothing_length` here is used for scaling - # TODO move this to the surface tension model since this is not a general thing cache.surface_normal[i, particle] += m_b / density_neighbor * - grad_kernel[i] * smoothing_length + grad_kernel[i] end cache.neighbor_count[particle] += 1 end @@ -112,8 +124,7 @@ function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system, return system end -function calc_normal_akinci!(system, neighbor_system, u_system, v, v_neighbor_system, - u_neighbor_system, semi, surfn) +function remove_invalid_normals!(system, surface_tension) # Normal not needed return system end @@ -132,18 +143,36 @@ function remove_invalid_normals!(system::FluidSystem, surface_tension::SurfaceTe return system end -function remove_invalid_normals!(system, surface_tension) - # Normal not needed +function remove_invalid_normals!(system::FluidSystem, surface_tension::SurfaceTensionMorris) + (; cache, smoothing_length) = system + + # We remove invalid normals i.e. they have a small norm (eq. 20) + normal_condition2 = (0.01 / smoothing_length)^2 + + for particle in each_moving_particle(system) + particle_surface_normal = cache.surface_normal[1:ndims(system), particle] + norm2 = dot(particle_surface_normal, particle_surface_normal) + + # see eq. 21 + if norm2 > normal_condition2 + cache.surface_normal[1:ndims(system), particle] = particle_surface_normal / + sqrt(norm2) + else + cache.surface_normal[1:ndims(system), particle] .= 0 + end + end + return system end -function compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) +function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) return system end -function compute_surface_normal!(system::FluidSystem, surface_tension::SurfaceTensionAkinci, +function compute_surface_normal!(system::FluidSystem, + surface_normal_method::ColorfieldSurfaceNormal, v, u, v_ode, u_ode, semi, t) - (; cache, surface_normal_method) = system + (; cache, surface_tension) = system # Reset surface normal set_zero!(cache.surface_normal) @@ -153,9 +182,69 @@ function compute_surface_normal!(system::FluidSystem, surface_tension::SurfaceTe u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) - calc_normal_akinci!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method) + calc_normal!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method) remove_invalid_normals!(system, surface_tension) end return system end + +# Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +function calc_curvature!(system::FluidSystem, neighbor_system::FluidSystem, u_system, v, + v_neighbor_system, u_neighbor_system, semi, surfn) + (; cache) = system + (; smoothing_kernel, smoothing_length) = surfn + + system_coords = current_coordinates(u_system, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + nhs = get_neighborhood_search(system, neighbor_system, semi) + + if smoothing_length != system.smoothing_length || + smoothing_kernel !== system.smoothing_kernel + # TODO: this is really slow but there is no way to easily implement multiple search radia + search_radius = compact_support(smoothing_kernel, smoothing_length) + nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius, + nparticles(system)) + PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords) + end + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, + nhs) do particle, neighbor, pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + density_neighbor = particle_density(v_neighbor_system, + neighbor_system, neighbor) + grad_kernel = kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) + for i in 1:ndims(system) + # The `smoothing_length` here is used for scaling + cache.surface_normal[i, particle] += m_b / density_neighbor * + grad_kernel[i] + end + + cache.neighbor_count[particle] += 1 + end + + return system +end + +function compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + return system +end + +function compute_curvature!(system::FluidSystem, + surface_tension::ViscosityMorris, + v, u, v_ode, u_ode, semi, t) + (; cache, surface_tension) = system + + # Reset surface curvature + set_zero!(cache.curvature) + + @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + + calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method) + end + return system +end diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 2d1f8e8f8..e3d991086 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -1,4 +1,5 @@ -abstract type AkinciTypeSurfaceTension end +abstract type SurfaceTensionModel end +abstract type AkinciTypeSurfaceTension <: SurfaceTensionModel end @doc raw""" CohesionForceAkinci(surface_tension_coefficient=1.0) @@ -50,6 +51,23 @@ struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension end end +struct SurfaceTensionMorris{ELTYPE} <: SurfaceTensionModel + surface_tension_coefficient::ELTYPE + + function SurfaceTensionMorris(; surface_tension_coefficient=1.0) + new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + end +end + +function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles) + return (;) +end + +function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) + curvature = Array{ELTYPE, 1}(undef, nparticles) + return (; curvature) +end + # Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`. # Currently, specializations reducing this to simple multiplications exist only up # to a power of three, see @@ -101,11 +119,18 @@ end return adhesion_force end +# Skip +@inline function surface_tension_force(surface_tension_a, surface_tension_b, + particle_system, neighbor_system, particle, neighbor, + pos_diff, distance, rho_a, rho_b) + return zero(pos_diff) +end + @inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, surface_tension_b::CohesionForceAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b) (; smoothing_length) = particle_system # No cohesion with oneself distance < sqrt(eps()) && return zero(pos_diff) @@ -120,7 +145,7 @@ end surface_tension_b::SurfaceTensionAkinci, particle_system::FluidSystem, neighbor_system::FluidSystem, particle, neighbor, - pos_diff, distance) + pos_diff, distance, rho_a, rho_b) (; smoothing_length, smoothing_kernel) = particle_system (; surface_tension_coefficient) = surface_tension_a @@ -134,14 +159,23 @@ end return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) .- - (surface_tension_coefficient * (n_a - n_b)) + (surface_tension_coefficient * (n_a - n_b) * smoothing_length) end -# Skip -@inline function surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, particle, neighbor, - pos_diff, distance) - return zero(pos_diff) +@inline function surface_tension_force(surface_tension_a::SurfaceTensionMorris, + surface_tension_b::SurfaceTensionMorris, + particle_system::FluidSystem, + neighbor_system::FluidSystem, particle, neighbor, + pos_diff, distance, rho_a, rho_b) + (; surface_tension_coefficient) = surface_tension_a + + # No surface tension with oneself + distance < sqrt(eps()) && return zero(pos_diff) + + n_a = surface_normal(particle_system, particle) + curvature_a = curvature(particle_system, particle) + + return surface_tension_coefficient / rho_a * curvature_a * n_a end @inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index bc0362a26..4a7ec4c5d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -61,7 +61,8 @@ function interact!(dv, v_particle_system, u_particle_system, dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, - particle, neighbor, pos_diff, distance) + particle, neighbor, pos_diff, distance, + rho_a, rho_b) dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 27a145115..499214927 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -119,7 +119,10 @@ function WeaklyCompressibleSPHSystem(initial_condition, create_cache_wcsph(correction, initial_condition.density, NDIMS, n_particles)..., cache...) cache = (; - create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., + create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, + n_particles)..., + create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, + n_particles)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, @@ -156,12 +159,6 @@ function create_cache_wcsph(::MixedKernelGradientCorrection, density, NDIMS, n_p return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end -function create_cache_wcsph(::SurfaceTensionAkinci, ELTYPE, NDIMS, nparticles) - surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) - neighbor_count = Array{ELTYPE, 1}(undef, nparticles) - return (; surface_normal, neighbor_count) -end - function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time @@ -240,7 +237,7 @@ function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_tension) = system + (; density_calculator, correction, surface_normal_method, surface_tension) = system compute_correction_values!(system, correction, u, v_ode, u_ode, semi) @@ -250,7 +247,8 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) compute_pressure!(system, v) - compute_surface_normal!(system, surface_tension, v, u, v_ode, u_ode, semi, t) + compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) return system end @@ -378,3 +376,8 @@ end @inline function surface_tension_model(system) return nothing end + +@inline function curvature(particle_system::FluidSystem, particle) + (; cache) = particle_system + return extract_svector(cache.curvature, particle_system, particle) +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index faee5ae6d..b8dab29a8 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -248,6 +248,10 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["neighbor_count"] = system.cache.neighbor_count end + if system.surface_tension isa ViscosityMorris + vtk["curvature"] = system.cache.curvature + end + if write_meta_data vtk["acceleration"] = system.acceleration vtk["viscosity"] = type2string(system.viscosity)