diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 29384b6d1..7beb184ed 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -55,13 +55,17 @@ of the boundary particle ``b``. ### Hydrodynamic density of dummy particles -We provide five options to compute the boundary density and pressure, determined by the `density_calculator`: +We provide six options to compute the boundary density and pressure, determined by the `density_calculator`: 1. (Recommended) With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the fluid according to [Adami et al., 2012](@cite Adami2012), and the density is obtained by applying the inverse of the state equation. This option usually yields the best results of the options listed here. -2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, +2. (Only relevant for FSI) With [`BernoulliPressureExtrapolation`](@ref), the pressure is extrapolated from the + pressure similar to the [`AdamiPressureExtrapolation`](@ref), but a relative velocity-dependent pressure part + is calculated between moving solids and fluids, which increases the boundary pressure in areas prone to + penetrations. +3. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, and the pressure is computed from the density with the state equation. -3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, +4. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, and the pressure is computed from the density with the state equation. Note that this causes a gap between fluid and boundary where the boundary is initialized without any contact to the fluid. This is due to overestimation of the boundary density @@ -69,10 +73,10 @@ We provide five options to compute the boundary density and pressure, determined contact to the fluid. Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat. See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13). -4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure +5. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure is computed from the density with the state equation. This option is not recommended. The other options yield significantly better results. -5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure +6. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure is not used. Instead, the fluid pressure is mirrored as boundary pressure in the momentum equation. This option is not recommended due to stability issues. See [`PressureMirroring`](@ref) @@ -93,7 +97,20 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den AdamiPressureExtrapolation ``` -#### 4. [`PressureZeroing`](@ref) +#### 2. [`BernoulliPressureExtrapolation`](@ref) +Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation: +```math +p_b = \frac{\sum_f (p_f + \frac{1}{2} \, \rho_{\text{neighbor}} \left( \frac{ (\mathbf{v}_f - \mathbf{v}_{\text{body}}) \cdot (\mathbf{x}_f - \mathbf{x}_{\text{neighbor}}) }{ \left\| \mathbf{x}_f - \mathbf{x}_{\text{neighbor}} \right\| } \right)^2 \times \text{factor} +\rho_f (\bm{g} - \bm{a}_b) \cdot \bm{r}_{bf}) W(\Vert r_{bf} \Vert, h)}{\sum_f W(\Vert r_{bf} \Vert, h)} +``` +where ``\mathbf{v}_f`` is the velocity of the fluid and ``\mathbf{v}_{\text{body}}`` is the velocity of the body. +This adjustment provides a higher boundary pressure for solid bodies moving with a relative velocity to the fluid to prevent penetration. +This modification is original and not derived from any literature source. + +```@docs + BernoulliPressureExtrapolation +``` + +#### 5. [`PressureZeroing`](@ref) This is the simplest way to implement dummy boundary particles. The density of each particle is set to the reference density and the pressure to the @@ -102,7 +119,7 @@ reference pressure (the corresponding pressure to the reference density by the s PressureZeroing ``` -#### 5. [`PressureMirroring`](@ref) +#### 6. [`PressureMirroring`](@ref) Instead of calculating density and pressure for each boundary particle, we modify the momentum equation, diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 6a2397a4b..a06056f8f 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -63,7 +63,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary -boundary_density_calculator = AdamiPressureExtrapolation() +boundary_density_calculator = BernoulliPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 89026040c..7f2d2272b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -70,7 +70,9 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka + PressureMirroring, PressureZeroing, BoundaryModelLastiwka, + BernoulliPressureExtrapolation + export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 1cd8e0dfe..8b7e507ed 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -70,11 +70,43 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, end @doc raw""" - AdamiPressureExtrapolation() + AdamiPressureExtrapolation(; pressure_offset=0.0) `density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. + """ -struct AdamiPressureExtrapolation end +struct AdamiPressureExtrapolation{ELTYPE} + pressure_offset::ELTYPE + + function AdamiPressureExtrapolation(; pressure_offset=0.0) + return new{eltype(pressure_offset)}(pressure_offset) + end +end + +@doc raw""" + BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=1.0) + +`density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. +- `factor=1.0` : Setting `factor` allows to just increase the strength of the dynamic + pressure part. + +""" +struct BernoulliPressureExtrapolation{ELTYPE} + pressure_offset :: ELTYPE + factor :: ELTYPE + + function BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=0.0) + return new{eltype(pressure_offset)}(pressure_offset, factor) + end +end @doc raw""" PressureMirroring() @@ -134,7 +166,9 @@ end @inline create_cache_model(initial_density, ::ContinuityDensity) = (; initial_density) -function create_cache_model(initial_density, ::AdamiPressureExtrapolation) +function create_cache_model(initial_density, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}) density = copy(initial_density) volume = similar(initial_density) @@ -194,7 +228,7 @@ end # Note that the other density calculators are dispatched in `density_calculators.jl` @inline function particle_density(v, ::Union{AdamiPressureExtrapolation, PressureMirroring, - PressureZeroing}, + PressureZeroing, BernoulliPressureExtrapolation}, boundary_model, particle) (; cache) = boundary_model @@ -216,10 +250,11 @@ end function compute_density!(boundary_model, ::Union{ContinuityDensity, AdamiPressureExtrapolation, + BernoulliPressureExtrapolation, PressureMirroring, PressureZeroing}, system, v, u, v_ode, u_ode, semi) # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`. - # For `AdamiPressureExtrapolation`, the density is updated in `compute_pressure!`. + # For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. return boundary_model end @@ -299,8 +334,10 @@ end boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) end -function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, - v_ode, u_ode, semi) +function compute_pressure!(boundary_model, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}, + system, v, u, v_ode, u_ode, semi) (; pressure, cache, viscosity) = boundary_model set_zero!(pressure) @@ -327,16 +364,18 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, nhs = get_neighborhood_search(neighbor_system, system, semi) # Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) else nhs = get_neighborhood_search(system, neighbor_system, semi) # Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) end @threaded system for particle in eachparticle(system) @@ -378,67 +417,81 @@ function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZe return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system, system_coords, - neighbor_coords, v_neighbor_system, - neighborhood_search) +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, system_coords, + neighbor_coords, v, + v_neighbor_system, + neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, - neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, + v, v_neighbor_system, + neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator foreach_point_neighbor(neighbor_system, system, neighbor_coords, system_coords, neighborhood_search; points=eachparticle(neighbor_system), parallel=false) do neighbor, particle, pos_diff, distance # Since neighbor and particle are switched pos_diff = -pos_diff - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) +@inline function boundary_pressure_extrapolation!(boundary_model, system, neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator # Loop over all pairs of particles and neighbors within the kernel cutoff. foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, neighborhood_search; points=eachparticle(system)) do particle, neighbor, pos_diff, distance - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) +@inline function boundary_pressure_inner!(boundary_model, boundary_density_calculator, + system, neighbor_system::FluidSystem, v, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure, + pressure_offset) density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) - resulting_acc = neighbor_system.acceleration - - current_acceleration(system, particle) + resulting_acceleration = neighbor_system.acceleration - + current_acceleration(system, particle) kernel_weight = smoothing_kernel(boundary_model, distance) - pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, - neighbor) + - dot(resulting_acc, density_neighbor * pos_diff)) * + pressure[particle] += (pressure_offset + + + particle_pressure(v_neighbor_system, neighbor_system, + neighbor) + + + dynamic_pressure(boundary_density_calculator, density_neighbor, + v, v_neighbor_system, pos_diff, distance, + particle, neighbor, system, neighbor_system) + + + dot(resulting_acceleration, density_neighbor * pos_diff)) * kernel_weight cache.volume[particle] += kernel_weight @@ -447,14 +500,46 @@ end kernel_weight, particle, neighbor) end +@inline function dynamic_pressure(boundary_density_calculator, density_neighbor, v, + v_neighbor_system, pos_diff, distance, particle, neighbor, + system, neighbor_system) + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::BoundarySystem, neighbor_system) + if system.ismoving[] + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) + + return 0.5 * boundary_density_calculator.factor * density_neighbor * + normal_velocity^2 / distance + end + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::SolidSystem, neighbor_system) + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) / distance + + return boundary_density_calculator.factor * density_neighbor * + dot(normal_velocity, normal_velocity) / 2 +end + function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) return cache end -function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, - neighbor_system, v_neighbor_system, kernel_weight, - particle, neighbor) +function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, neighbor_system, + v_neighbor_system, kernel_weight, particle, neighbor) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) for dim in 1:ndims(neighbor_system) diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 27caa0293..a72ac409c 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -54,12 +54,14 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, - v_fluid .* - ones(size(fluid.coordinates)), - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid .* + ones(size(fluid.coordinates)), + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps() @@ -92,10 +94,12 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, v_fluid, - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid, + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps()