From 0f3536c8074e46a6436301df3f1cedc7aa7a5a6f Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 15 Jul 2024 17:29:48 +0200 Subject: [PATCH] add new bnd density calculator --- docs/src/systems/boundary.md | 14 +- examples/fsi/falling_spheres_2d.jl | 2 +- src/TrixiParticles.jl | 2 +- src/general/system.jl | 5 + .../dummy_particles/dummy_particles.jl | 129 ++++++++++++------ 5 files changed, 106 insertions(+), 46 deletions(-) diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index ae4b9b676..f2cf29d64 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), 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`, the pressure is extrapolated from the + pressure similar to the `AdamiPressureExtrapolation`, 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) 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 846bfa4a0..743fb74e1 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,7 +65,7 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing + PressureMirroring, PressureZeroing, BernoulliPressureExtrapolation export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/general/system.jl b/src/general/system.jl index 5d3143043..49e95fbc8 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -93,6 +93,11 @@ end @inline current_velocity(v, system, particle) = extract_svector(v, system, particle) +@inline function current_acceleration(system, particle) + # TODO: Return `dv` of solid particles + return zero(SVector{ndims(system), eltype(system)}) +end + @inline function smoothing_kernel(system, distance) (; smoothing_kernel, smoothing_length) = system return kernel(smoothing_kernel, distance, smoothing_length) diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index c056e465c..2b132de99 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -87,6 +87,24 @@ struct AdamiPressureExtrapolation{ELTYPE} end end +@doc raw""" + BernoulliPressureExtrapolation(; 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 BernoulliPressureExtrapolation{ELTYPE} + pressure_offset::ELTYPE + + function BernoulliPressureExtrapolation(; pressure_offset=0.0) + return new{eltype(pressure_offset)}(pressure_offset) + end +end + @doc raw""" PressureMirroring() @@ -145,7 +163,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) @@ -205,7 +225,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 @@ -227,10 +247,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 @@ -310,7 +331,10 @@ end boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) end -function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, +function compute_pressure!(boundary_model, + bnd_density_calc::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}, system, + v, u, v_ode, u_ode, semi) (; pressure, cache, viscosity) = boundary_model @@ -338,16 +362,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, - v_neighbor_system, nhs) + bnd_pressure_extrapolation_neighbor!(boundary_model, bnd_density_calc, 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, - v_neighbor_system, nhs) + bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) end @threaded system for particle in eachparticle(system) @@ -406,43 +432,43 @@ end v_neighbor_system[1:ndims(system), neighbor])^2 end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system, system_coords, - neighbor_coords, v_neighbor_system, - neighborhood_search) +@inline function bnd_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, system_coords, + neighbor_coords, 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, - v_neighbor_system, - neighborhood_search) +@inline function bnd_pressure_extrapolation_neighbor!(boundary_model, bnd_density_calc, + system, neighbor_system::FluidSystem, + system_coords, neighbor_coords, v, + v_neighbor_system, + neighborhood_search) (; pressure, cache, viscosity) = 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, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure, pressure_offset) + bnd_pressure_inner!(boundary_model, bnd_density_calc, 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, - v_neighbor_system, neighborhood_search) +@inline function bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, 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, - v_neighbor_system, neighborhood_search) +@inline function bnd_pressure_extrapolation!(boundary_model, bnd_density_calc, 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 @@ -451,19 +477,21 @@ end neighborhood_search; points=eachparticle(system)) do particle, neighbor, pos_diff, distance - adami_pressure_inner!(boundary_model, system, neighbor_system, v, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure, pressure_offset) + bnd_pressure_inner!(boundary_model, bnd_density_calc, 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::FluidSystem, v, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure, - pressure_offset) +@inline function bnd_pressure_inner!(boundary_model, ::BernoulliPressureExtrapolation, + 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) + kernel_weight = smoothing_kernel(boundary_model, distance) pressure[particle] += (pressure_offset + @@ -472,7 +500,30 @@ end dynamic_pressure(density_neighbor, v, v_neighbor_system, particle, neighbor, system) + - dot(neighbor_system.acceleration, + dot(resulting_acc, + density_neighbor * pos_diff)) * kernel_weight + + cache.volume[particle] += kernel_weight + + compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, + kernel_weight, particle, neighbor) +end + +@inline function bnd_pressure_inner!(boundary_model, ::AdamiPressureExtrapolation, + 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) + + kernel_weight = smoothing_kernel(boundary_model, distance) + + pressure[particle] += (pressure_offset + + particle_pressure(v_neighbor_system, neighbor_system, + neighbor) + + dot(resulting_acc, density_neighbor * pos_diff)) * kernel_weight cache.volume[particle] += kernel_weight