diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 59f93c0766..0cd90d517d 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -220,9 +220,9 @@ Pages = [joinpath("general", "corrections.jl")] ### Overview of surface normal calculation in SPH -Surface normals are essential for modeling surface tension as they provide the directionality -of forces acting at the fluid interface. They are calculated based on the particle properties and -their spatial distribution. +Surface normals provide the directionality of forces acting at the fluid interface. They are +used by the full Akinci model and both Morris models, but not by the cohesion-only Akinci model. +They are calculated based on the particle properties and their spatial distribution. #### Color field and gradient-based surface normals @@ -296,6 +296,32 @@ In the following table some values are shown for reference. The values marked wi | **Water** | 0.07288 [Lange](@cite Lange2005) | | **Mercury** | 0.486502 [Lange](@cite Lange2005) | +### Model configuration + +All surface tension coefficients must be finite and non-negative. A zero coefficient disables +the fluid-fluid surface force. Wall adhesion is controlled independently by the boundary's +`adhesion_coefficient`. + +`CohesionForceAkinci` only evaluates the pairwise cohesion and optional wall-adhesion forces. +It does not require surface normals or `reference_particle_spacing`. The full +`SurfaceTensionAkinci` model and both Morris models require a surface-normal method. When one +of these models is selected without an explicit method, `ColorfieldSurfaceNormal()` is used. +Runnable configurations are available in `examples/fluid/cohesion_force_akinci_2d.jl` and +`examples/fluid/akinci_wetting_2d.jl`. + +!!! warning "Akinci coefficients in two dimensions" + The Akinci cohesion and adhesion kernels implemented here use the normalization published + for the three-dimensional model. In two-dimensional simulations, their coefficients are + empirical numerical parameters rather than resolution-independent physical values in N/m. + Recheck the coefficient when changing particle spacing or smoothing length. + + Both published kernels have dimensions of inverse volume: the cohesion numerator scales + with ``L^6`` and its denominator with ``L^9``, while the fourth-root term in the adhesion + kernel changes its ``L^{-3.25}`` prefactor to ``L^{-3}``. Particle mass scales with + ``L^D`` in ``D`` dimensions. At a fixed smoothing-length-to-spacing ratio, the resulting + pair contribution therefore scales with ``L^{D-3}``: it is resolution-independent in 3D + but retains an inverse-length factor in 2D. + ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf) The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components, diff --git a/examples/fluid/akinci_wetting_2d.jl b/examples/fluid/akinci_wetting_2d.jl new file mode 100644 index 0000000000..e15c5c1c8c --- /dev/null +++ b/examples/fluid/akinci_wetting_2d.jl @@ -0,0 +1,67 @@ +# ========================================================================================== +# 2D Akinci Surface Tension and Wall Adhesion +# +# A circular drop rests on the bottom wall. Set `wetting=true` to compare stronger wall +# adhesion with a surface-tension-dominated non-wetting setup. The full Akinci model uses +# surface normals, while its wall interaction uses the Akinci adhesion kernel. +# ========================================================================================== + +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +particle_spacing = 0.005 +fluid_density = 1000.0 +sound_speed = 100.0 +gravity = 9.81 +tspan = (0.0, 0.3) +wetting = false + +if wetting + surface_tension_coefficient = 0.01 + adhesion_coefficient = 1.0 + nu = 0.0005 +else + surface_tension_coefficient = 2.0 + adhesion_coefficient = 0.001 + nu = 0.001 +end + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1) + +tank = RectangularTank(particle_spacing, (0.0, 0.0), (0.5, 0.1), fluid_density; + n_layers=3, faces=(true, true, true, false), + acceleration=(0.0, -gravity), state_equation) +drop = SphereShape(particle_spacing, 0.05, (0.25, 0.05), fluid_density; + sphere_type=VoxelSphere()) + +smoothing_length = particle_spacing - eps() +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +alpha = 8 * nu / (smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +fluid_system = WeaklyCompressibleSPHSystem(drop; smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), + state_equation, viscosity, + acceleration=(0.0, -gravity), + surface_tension=SurfaceTensionAkinci(; + surface_tension_coefficient), + correction=AkinciFreeSurfaceCorrection(fluid_density), + reference_particle_spacing=particle_spacing) + +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system, + boundary_density_calculator=AdamiPressureExtrapolation(), + viscosity=ViscosityAdami(; nu=4 * nu), + clip_negative_pressure=true) +boundary_system = WallBoundarySystem(tank.boundary, boundary_model; + adhesion_coefficient) + +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.01) +callbacks = CallbackSet(info_callback, saving_callback) + +sol = solve(ode, RDPK3SpFSAL35(); abstol=1e-7, reltol=1e-4, + save_everystep=false, callback=callbacks) diff --git a/examples/fluid/cohesion_force_akinci_2d.jl b/examples/fluid/cohesion_force_akinci_2d.jl new file mode 100644 index 0000000000..38527dbdf9 --- /dev/null +++ b/examples/fluid/cohesion_force_akinci_2d.jl @@ -0,0 +1,56 @@ +# ========================================================================================== +# 2D Cohesion-Only Akinci Surface Force +# +# This example evolves a rectangular fluid patch with `CohesionForceAkinci`. The +# cohesion-only model does not calculate surface normals and therefore does not require +# `reference_particle_spacing`. In 2D, its coefficient is an empirical numerical parameter. +# ========================================================================================== + +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +particle_spacing = 0.025 +fluid_size = (0.2, 0.1) +fluid_density = 1000.0 +sound_speed = 20.0 +tspan = (0.0, 0.2) + +fluid = RectangularShape(particle_spacing, + round.(Int, fluid_size ./ particle_spacing), + zeros(length(fluid_size)); density=fluid_density) + +smoothing_length = particle_spacing - eps() +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, clip_negative_pressure=true) + +nu = 0.01 +alpha = 8 * nu / (smoothing_length * sound_speed) +viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) +surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.001) + +fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=SummationDensity(), + state_equation, viscosity, surface_tension, + source_terms=SourceTermDamping(damping_coefficient=0.5)) + +semi = Semidiscretization(fluid_system) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02) +stepsize_callback = StepsizeCallback(cfl=0.5) +callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, save_everystep=false, callback=callbacks) + +v_ode, u_ode = sol.u[end].x +v = TrixiParticles.wrap_v(v_ode, fluid_system, semi) +velocity = TrixiParticles.current_velocity(v, fluid_system) +total_momentum = vec(sum(velocity .* transpose(fluid_system.mass); dims=2)) +center_of_mass_velocity = total_momentum / sum(fluid_system.mass) +final_kinetic_energy = kinetic_energy(fluid_system, nothing, nothing, v_ode, u_ode, semi, + sol.t[end]) + +@info "Cohesion diagnostics" center_of_mass_velocity final_kinetic_energy diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 22969a427c..8f6a387221 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -65,12 +65,15 @@ viscosity_fluid = ArtificialViscosityMonaghan(; alpha, beta=0.0) density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) # density_diffusion = DensityDiffusionAntuono(delta=0.1) +surface_tension = nothing +reference_particle_spacing = 0 + fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothing_length, density_calculator=fluid_density_calculator, state_equation, viscosity=viscosity_fluid, density_diffusion, acceleration=(0.0, -gravity), - correction=nothing, surface_tension=nothing, - reference_particle_spacing=0) + correction=nothing, surface_tension, + reference_particle_spacing) # ========================================================================================== # ==== Boundary diff --git a/examples/fluid/dam_break_oil_film_2d.jl b/examples/fluid/dam_break_oil_film_2d.jl index 7a890eb918..3d57557097 100644 --- a/examples/fluid/dam_break_oil_film_2d.jl +++ b/examples/fluid/dam_break_oil_film_2d.jl @@ -33,12 +33,15 @@ nu_sim_water = nu_ratio * nu_sim_oil oil_viscosity = ViscosityMorris(nu=nu_sim_oil) -# TODO: broken if both systems use surface tension +# A physically consistent two-phase surface-tension interaction requires an explicit +# interface model. Until that model is available, this example focuses on density and +# viscosity contrast and keeps surface tension disabled on both fluids. +surface_tension = nothing trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"); sol=nothing, fluid_particle_spacing, tspan, viscosity_fluid=ViscosityMorris(nu=nu_sim_water), smoothing_length, gravity, density_diffusion=nothing, sound_speed, prefix="", - reference_particle_spacing=fluid_particle_spacing) + surface_tension) # ========================================================================================== # ==== Setup oil layer @@ -65,18 +68,7 @@ oil_system = WeaklyCompressibleSPHSystem(oil; state_equation=oil_eos, viscosity=oil_viscosity, acceleration=(0.0, -gravity), - surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01), - correction=AkinciFreeSurfaceCorrection(oil_density), - reference_particle_spacing=fluid_particle_spacing) - -# oil_system = WeaklyCompressibleSPHSystem(oil; -# smoothing_kernel, smoothing_length, -# density_calculator=fluid_density_calculator, -# state_equation=oil_eos, -# viscosity=oil_viscosity, -# acceleration=(0.0, -gravity), -# surface_tension=SurfaceTensionMorris(surface_tension_coefficient=0.03), -# reference_particle_spacing=fluid_particle_spacing) + surface_tension) # ========================================================================================== # ==== Simulation diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 657f496a48..ccdc876865 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -1,9 +1,8 @@ # ========================================================================================== -# 2D Falling Water Spheres Simulation (With and Without Surface Tension) +# 2D Falling Water Spheres Simulation with Surface Tension # # This example simulates two circular water "spheres" falling under gravity. -# One sphere includes a surface tension model (Akinci et al.), while the other does not. -# This demonstrates the effect of surface tension on fluid behavior. +# Both spheres belong to one fluid system and use the same surface tension model. # ========================================================================================== using TrixiParticles @@ -44,6 +43,7 @@ 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)) +falling_spheres = isnothing(sphere2) ? sphere1 : union(sphere1, sphere2) # ========================================================================================== # ==== Fluid @@ -53,28 +53,20 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center, 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, beta=0.0) -density_diffusion = DensityDiffusionAntuono(delta=0.1) surface_tension_coefficient = 0.05 surface_tension = SurfaceTensionAkinci(; surface_tension_coefficient) +reference_particle_spacing = isnothing(surface_tension) ? 0 : fluid_particle_spacing -sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1; +sphere_surface_tension = EntropicallyDampedSPHSystem(falling_spheres; smoothing_kernel=fluid_smoothing_kernel, smoothing_length=fluid_smoothing_length, sound_speed, viscosity, density_calculator=ContinuityDensity(), acceleration, surface_tension, - reference_particle_spacing=fluid_particle_spacing) - -sphere = WeaklyCompressibleSPHSystem(sphere2; smoothing_kernel=fluid_smoothing_kernel, - smoothing_length=fluid_smoothing_length, - density_calculator=fluid_density_calculator, - state_equation, viscosity, density_diffusion, - acceleration) + reference_particle_spacing) # ========================================================================================== # ==== Boundary @@ -94,7 +86,7 @@ boundary_system = WallBoundarySystem(tank.boundary, boundary_model; # ========================================================================================== # ==== Simulation -semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system) +semi = Semidiscretization(sphere_surface_tension, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=1000) diff --git a/examples/fluid/falling_water_spheres_3d.jl b/examples/fluid/falling_water_spheres_3d.jl index 511927145f..84db78fe13 100644 --- a/examples/fluid/falling_water_spheres_3d.jl +++ b/examples/fluid/falling_water_spheres_3d.jl @@ -1,10 +1,8 @@ # ========================================================================================== -# 3D Falling Water Spheres Simulation (With and Without Surface Tension) +# 3D Falling Water Spheres Simulation with Surface Tension # # This example extends `falling_water_spheres_2d.jl` to three dimensions. -# It simulates two spherical volumes of water falling under gravity. -# One sphere includes a surface tension model, while the other does not, -# demonstrating the effect of surface tension in 3D. +# It simulates two spherical volumes of water in one fluid system falling under gravity. # ========================================================================================== using TrixiParticles diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl index 1abcdcb149..36baaeab07 100644 --- a/examples/fluid/sphere_surface_tension_2d.jl +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -23,7 +23,7 @@ sound_speed = 20.0 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7, clip_negative_pressure=true) -# For all surface tension simulations, we need a compact support of `2 * particle_spacing` +# The surface tension configurations below use a compact support of `2 * particle_spacing`. # smoothing_length = particle_spacing # smoothing_kernel = WendlandC2Kernel{2}() # nu = 0.01 diff --git a/examples/fluid/sphere_surface_tension_3d.jl b/examples/fluid/sphere_surface_tension_3d.jl index 7cc8ce78af..6e759bc3a3 100644 --- a/examples/fluid/sphere_surface_tension_3d.jl +++ b/examples/fluid/sphere_surface_tension_3d.jl @@ -17,7 +17,7 @@ fluid_size = (0.9, 0.9, 0.9) sound_speed = 20.0 -# For all surface tension simulations, we need a compact support of `2 * particle_spacing` +# The surface tension configurations below use a compact support of `2 * particle_spacing`. smoothing_length = 1.0 * particle_spacing nu = 0.04 diff --git a/examples/fluid/sphere_surface_tension_wall_2d.jl b/examples/fluid/sphere_surface_tension_wall_2d.jl index 5065536c93..b651601de9 100644 --- a/examples/fluid/sphere_surface_tension_wall_2d.jl +++ b/examples/fluid/sphere_surface_tension_wall_2d.jl @@ -63,6 +63,6 @@ sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1; reference_particle_spacing=fluid_particle_spacing) trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"); - sphere=nothing, sphere1, adhesion_coefficient=0.001, wall_viscosity=4.0 * nu, + sphere1, sphere2=nothing, adhesion_coefficient=0.001, wall_viscosity=4.0 * nu, alpha, sound_speed, fluid_density, nu, fluid_particle_spacing, tspan, tank_size, fluid_smoothing_length, sphere_surface_tension) diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index e0ec835c1e..75ab98624e 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -33,7 +33,7 @@ Boundary model for [`WallBoundarySystem`](@ref). in areas of low pressure, against which the particle shifting technique is fighting. - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. + which is needed when using a surface-normal method. # Examples ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1) # Free-slip condition diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 3720b57ba1..baa2215b0f 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -51,9 +51,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. - `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) - `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) + (default: no surface normal method or `ColorfieldSurfaceNormal()` + if the surface tension model requires normals) - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. + which is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation @@ -119,12 +120,11 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end + surface_normal_method = default_surface_normal_method(surface_tension, + surface_normal_method) if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method")) end if correction isa ShepardKernelCorrection && diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 777480c923..0d42521148 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -229,10 +229,13 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, sem if surface_tension isa SurfaceTensionMorris || surface_tension isa SurfaceTensionMomentumMorris - v = wrap_v(v_ode, system, semi) - dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 / - (2 * pi * surface_tension.surface_tension_coefficient)) - dt = min(dt, dt_surface_tension) + coefficient = surface_tension.surface_tension_coefficient + if !iszero(coefficient) + v = wrap_v(v_ode, system, semi) + dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 / + (2 * pi * coefficient)) + dt = min(dt, dt_surface_tension) + end end return dt diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 4db94ea763..8a7a72e47f 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -17,8 +17,17 @@ end function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, ideal_density_threshold=0.0) - return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, - ideal_density_threshold) + thresholds = promote(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + return ColorfieldSurfaceNormal(thresholds...) +end + +@inline function default_surface_normal_method(surface_tension, surface_normal_method) + if isnothing(surface_normal_method) && requires_surface_normal(surface_tension) + return ColorfieldSurfaceNormal() + end + + return surface_normal_method end function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..04067e5087 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -1,22 +1,38 @@ abstract type AbstractSurfaceTension end abstract type AkinciTypeSurfaceTension <: AbstractSurfaceTension end +function validate_surface_tension_coefficient(surface_tension_coefficient) + if !(surface_tension_coefficient isa Real) || + !isfinite(surface_tension_coefficient) || surface_tension_coefficient < 0 + throw(ArgumentError("`surface_tension_coefficient` must be a finite, non-negative real number")) + end + + return surface_tension_coefficient +end + @doc raw""" CohesionForceAkinci(surface_tension_coefficient=1.0) This model only implements the cohesion force of the Akinci [Akinci2013](@cite) surface tension model. +It does not require a surface-normal method. + +The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional +simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and +may need to be adjusted when changing the resolution. See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: Modifies the intensity of the surface tension-induced force, - enabling the tuning of the fluid's surface tension properties within the simulation. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient modifying the + fluid-fluid cohesion force. Zero disables this force; wall adhesion is controlled by the + boundary's `adhesion_coefficient`. """ -struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension +struct CohesionForceAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function CohesionForceAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -28,18 +44,22 @@ principles outlined by Akinci [Akinci2013](@cite). This model is instrumental in behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or separation, by utilizing intra-particle forces. +The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional +simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and +may need to be adjusted when changing the resolution. + See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the magnitude of - surface tension forces, facilitating the fine-tuning of how surface tension phenomena - are represented in the simulation. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + magnitude of surface tension forces. Zero disables the fluid-fluid force. """ -struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension +struct SurfaceTensionAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -55,14 +75,15 @@ See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: Adjusts the magnitude of the surface tension - forces, enabling tuning of fluid surface behaviors in simulations. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + magnitude of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMorris{ELTYPE} <: AbstractSurfaceTension +struct SurfaceTensionMorris{ELTYPE <: Real} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -87,17 +108,24 @@ numerical adjustments at higher resolutions. See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the strength of surface tension - forces, allowing fine-tuning to replicate physical behavior. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + strength of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension +struct SurfaceTensionMomentumMorris{ELTYPE <: Real} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end +# Surface-model capabilities are expressed through dispatch so that constructors and update +# stages do not need to duplicate concrete model checks. +@inline requires_surface_normal(::Nothing) = false +@inline requires_surface_normal(::CohesionForceAkinci) = false +@inline requires_surface_normal(::Any) = true + function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) @@ -122,6 +150,8 @@ end (; surface_tension_coefficient) = surface_tension # Eq. 2 + # This is the three-dimensional normalization published by Akinci et al. In 2D, the + # coefficient is an empirical numerical parameter; see the model docstring. # We only reach this function when `sqrt(eps()) < distance <= support_radius` if distance > 0.5 * support_radius # Attractive force @@ -141,18 +171,17 @@ end @inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, adhesion_coefficient) - - # The neighborhood search has an `<=` check, but for `distance == support_radius` - # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`. - # TODO Change this in the neighborhood search? - # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19 distance >= support_radius && return zero(pos_diff) distance <= 0.5 * support_radius && return zero(pos_diff) - # Eq. 7 - A = 0.007 / support_radius^3.25 * - (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25 + # Eq. 7. The factored radicand avoids cancellation close to the support boundary. + radicand = 2 * (2 * distance - support_radius) * + (support_radius - distance) / support_radius + fourth_root = sqrt(sqrt(max(zero(radicand), radicand))) + normalization = convert(typeof(support_radius), 0.007) / + (support_radius^3 * sqrt(sqrt(support_radius))) + A = normalization * fourth_root # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as # `m_b = rho_0 * volume` (Akinci boundary condition treatment) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eea0607d7d..c97c829824 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -54,9 +54,10 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. gravity-like source terms. - `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) - `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) + (default: no surface normal method or `ColorfieldSurfaceNormal()` + if the surface tension model requires normals) - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. + which is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation @@ -130,12 +131,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end + surface_normal_method = default_surface_normal_method(surface_tension, + surface_normal_method) if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method")) end pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index f7fc1890d6..1007df1f43 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -570,6 +570,37 @@ @test count_rhs_allocations(sol) == 0 end + @trixi_testset "fluid/cohesion_force_akinci_2d.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "cohesion_force_akinci_2d.jl"), + particle_spacing=0.05, tspan=(0.0, 0.01), + saving_callback=nothing) + @test sol.retcode == ReturnCode.Success + @test isnothing(fluid_system.surface_normal_method) + @test !haskey(fluid_system.cache, :surface_normal) + @test norm(center_of_mass_velocity) < 100eps() + @test final_kinetic_energy >= 0 + @test count_rhs_allocations(sol) == 0 + end + + @trixi_testset "fluid/akinci_wetting_2d.jl" begin + for wetting in (false, true) + @testset "wetting=$wetting" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "akinci_wetting_2d.jl"), + particle_spacing=0.02, + tspan=(0.0, 0.001), wetting, + saving_callback=nothing) + @test sol.retcode == ReturnCode.Success + @test fluid_system.surface_normal_method isa ColorfieldSurfaceNormal + @test adhesion_coefficient == (wetting ? 1.0 : 0.001) + @test count_rhs_allocations(sol) == 0 + end + end + end + @trixi_testset "fluid/sphere_surface_tension_2d.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 7fe8abbd97..d470da0ba9 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -1,5 +1,119 @@ - @testset verbose=true "Surface Tension" begin + @testset "constructors and capabilities" begin + constructors = (CohesionForceAkinci, SurfaceTensionAkinci, + SurfaceTensionMorris, SurfaceTensionMomentumMorris) + + for constructor in constructors + model = constructor(surface_tension_coefficient=0.5f0) + @test model.surface_tension_coefficient === 0.5f0 + @test iszero(constructor(surface_tension_coefficient=0).surface_tension_coefficient) + + for coefficient in (-1.0, NaN, Inf, -Inf, 1.0im, "invalid") + @test_throws ArgumentError constructor(surface_tension_coefficient=coefficient) + end + end + + @test !TrixiParticles.requires_surface_normal(nothing) + @test !TrixiParticles.requires_surface_normal(CohesionForceAkinci()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionAkinci()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionMorris()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionMomentumMorris()) + + normal_method = ColorfieldSurfaceNormal(boundary_contact_threshold=1, + interface_threshold=0.1f0, + ideal_density_threshold=0.25) + @test normal_method isa ColorfieldSurfaceNormal{Float64} + @test ColorfieldSurfaceNormal(boundary_contact_threshold=0.1f0, + interface_threshold=0.01f0, + ideal_density_threshold=0.0f0) isa + ColorfieldSurfaceNormal{Float32} + end + + @testset "cohesion-only systems do not require normals" begin + coordinates = [0.0 1.0; + 0.0 0.0] + initial_condition = InitialCondition(; coordinates, density=ones(2), + particle_spacing=1.0) + smoothing_kernel = WendlandC2Kernel{2}() + smoothing_length = 1.0 + surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.1) + + wcsph = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension) + edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension) + + for system in (wcsph, edac) + @test isnothing(system.surface_normal_method) + @test !haskey(system.cache, :surface_normal) + @test !haskey(system.cache, :neighbor_count) + @test !haskey(system.cache, :reference_particle_spacing) + end + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci()) + @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension=SurfaceTensionAkinci()) + + full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) + @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal + @test haskey(full_akinci.cache, :surface_normal) + end + + @testset "zero Morris coefficient does not restrict the time step" begin + function calculate_initial_dt(surface_tension) + initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], + density=ones(2), particle_spacing=1.0) + reference_particle_spacing = isnothing(surface_tension) ? 0 : 1.0 + system = WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel=WendlandC2Kernel{2}(), + smoothing_length=1.0, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension, + reference_particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.1)) + v_ode, u_ode = ode.u0.x + return TrixiParticles.calculate_dt(v_ode, u_ode, 0.25, semi.systems[1], semi) + end + + dt_without_surface_tension = calculate_initial_dt(nothing) + dt_with_zero_csf = calculate_initial_dt(SurfaceTensionMorris(; + surface_tension_coefficient=0.0)) + dt_with_zero_css = calculate_initial_dt(SurfaceTensionMomentumMorris(; + surface_tension_coefficient=0.0)) + + @test dt_with_zero_csf == dt_without_surface_tension + @test dt_with_zero_css == dt_without_surface_tension + end + @testset verbose=true "`cohesion_force_akinci`" begin surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) support_radius = 1.0 @@ -88,6 +202,16 @@ test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) + + support_radius_f32 = 15.594092f0 + distance_f32 = prevfloat(support_radius_f32) + near_support = TrixiParticles.adhesion_force_akinci(surface_tension, + support_radius_f32, 1.0f0, + Float32[1, 0], distance_f32, + 1.0f0) + @test eltype(near_support) == Float32 + @test all(isfinite, near_support) + @test 0 < norm(near_support) < eps(Float32) end @testset "compute_stress_tensors! (MomentumMorris)" begin