Decoupled Boundary Handling in SPH

Particle-based boundary representations are frequently used in Smoothed Particle Hydrodynamics (SPH) due to their simple integration into fluid solvers. Commonly, incompressible fluid solvers estimate the current density and corresponding forces in case the current density exceeds the rest density to push fluid particles apart. Close to the boundary, the calculation of the fluid particles' density involves both, neighboring fluid and neighboring boundary particles, yielding an overestimation of density, and, subsequently, wrong pressure forces and wrong velocities leading to the disturbed fluid particles' behavior in the vicinity of the boundary. In this paper, we present a detailed explanation of this disturbed fluid particle behavior, which is mainly due to the combined or coupled handling of the fluid-fluid particle and the fluid-boundary particle interaction. We propose the decoupled handling of both interaction types, leading to two densities for a given fluid particle, i.e., fluid-induced density and boundary-induced density. In our approach, we alternately apply the corresponding fluid-induced and boundary-induced forces during pressure estimation. This separation avoids force overestimation and reduces unintended fluid dynamics near the boundary, as well as a consistent fluid-boundary distance across different fluid amounts and different particle-based boundary handling methods. We compare our method with two regular state-of-the-art methods in different experiments and show how our method handles detailed boundary shapes.


Introduction
Boundary handling is an important part of fluid simulation using Smoothed Particle Hydrodynamics (SPH).There are several approaches to handle the interaction of SPH-based fluids with boundaries, i.e., particle-based Akinci et al. (2012), external boundary handling methods Losasso et al. (2008) or direct boundary integral methods Koschier and Bender (2017).
In particle-based boundary handling, a particle's density is usually calculated by taking all neighboring particles into account, i.e., fluid and boundary particles.Whereas including fluid particles in the neighborhood is a regular rule in SPH, including static boundary particles causes the wrong estimation of forces between fluid particles in the vicinity of boundaries, as we demonstrate in this paper.This yields distances between the fluid and the boundary that depend on the number of fluid particles and result in inconsistent fluid behavior and unintended fluid dynamics close to boundaries (see, for instance, Fig. 4).
In this paper, we first analyze the cause of the inconsistent fluid behavior close to boundaries, and we propose a new method for decoupled boundary density calculation, which divides the density calculation into the boundary-induced density and the fluid-induced density.While the first one takes only the neighboring boundary particles of a given fluid particle into account, the latter accounts for the interaction of a given fluid particle with its neighboring fluid particles.
We also propose a decoupled pressure calculation method, which alternates between the boundaryinduced density and the fluid-induced density.This separation drastically reduces the overestimation of forces appearing in the standard coupled version and makes the fluid quantities more consistent across the complete fluid domain, e.g., it reduces residual oscillations caused by force overestimation.
In our experiments, we integrate our method with the IISPH solver Ihmsen et al. (2013) (in the variant by Band et al. 2018a) using two types of boundary pressure estimators, i.e., Pressure Mirroring (PM) Akinci et al. (2012) and Pressure Boundaries (PB) Band et al. (2018a).Comparing the decoupled with the mixed methods, the decoupled approach demonstrates significant benefits in terms of residual oscillations, residual pressure, and velocity, moreover, the decoupled method also brings more consistency across different pressure estimators.
Conceptually, our decoupled approach can be combined with any scheme for estimating the density of the boundary particles, and it can be applied to nonparticle-based boundary handling schemes as well.

Prior Work
One of the first SPH approaches in fluid-boundary interaction is described by Monaghan Monaghan (1994), where the boundary is represented by particles and the interaction between fluid and boundary particles is ruled by the Lennard-Jones potential.Later, alternative methods of calculation of the forces between fluid and boundary particles have been proposed, such as direct forcing by Becker et al. Becker et al. (2009) or pressure forces by Akinci et al. 2012.Modern approaches in boundary handling estimate pressure values for the boundary particles that yield forces and resulting accelerations for fluid particles that prevent penetration of fluid particles into the boundary Koschier et al. (2020Koschier et al. ( , 2022)).
Several approaches to pressure calculation for boundary particles are used in computer graphics.The Pressure Mirroring (PM) method Akinci et al. (2012) estimates the pressure for a boundary particle by directly assigning the pressure of the currently interacting fluid particle.Other methods use a pressure extrapolation approach, where the boundary particles' pressure is extrapolated based on the Arrows represent an interaction between the selected fluid particle and its neighboring fluid particles (blue) and boundary particles (yellow).Blue lines represent interaction with fluid neighbors, and yellow arrows represent interaction with boundary particles.The left flowchart is the classical approach in density and pressure calculation, and the right flowchart depicts our decoupled approach pressure and position of the neighboring fluid particles.Specifically, the Pressure Boundary (PB) method uses a separate pressure solver to estimate a unique pressure value for each boundary particle Band et al. (2018a).An alternative approach uses a Moving Least Squares pressure extrapolation to calculate the pressure at boundary particles, where the pressure values of the boundary particles correspond to the pressure values of the fluid particles in its vicinity Band et al. (2018b).
Beyond particle-based boundary handling methods, other approaches have been proposed, such as direct forcing without boundary particles applied to fluid particles that approach the boundary surface, as done by Bodin et al. 2011.Other approaches to nonparticle-based boundary representation include the meshed-based method by Huber et al. 2015, the distance function-based method by Harada et al. 2007, the density maps by Bender et al. 2019, and the semianalytical boundary handling by Winchenbach et al. 2020.
Almost all of these boundary handling methods perform density calculations to deduce pressure forces that prevent boundary penetration.These methods estimate the density using both, boundaryand fluid-induced densities simultaneously in a coupled way, leading to an overestimation of pressure and, hence, forces, and leads to the excess motion of the fluid particles; see detailed discussion in Sec.3.
Our approach of decoupled density estimation borrows from a similar concept proposed by Gissler et al. 2019 applied in the context of boundary-boundary interactions, i.e., for rigid bodies collision reaction.For the calculation of regular SPH quantities, Gissler et al. 2019 use standard "coupled" density estimation, while they exclude fluid particles when estimating forces between rigid objects.

Problem Statement
The main challenge regarding boundary handling in SPH is the calculation of forces induced by the boundary on the fluid particles.The boundary must have at least the property of impenetrability, i.e., the forces induced by the boundary must be large enough to push the approaching fluid particles back to prevent penetration.
Fig. 1 depicts the interaction of a single fluid particle (highlighted with a dark blue outline) with neighboring fluid (blue) and boundary (yellow) particles, where blue and orange arrows represent interactions with fluid and boundary neighbors, respectively.Both types of interaction serve different goals, that is, the interaction with the boundary particles in the fluid particle's neighborhood provides impenetrability, while the interaction with the neighboring fluid particles enables fluid dynamics.
Commonly in SPH, the calculation of the density and the pressure forces for the given fluid particle considers all neighboring fluid and boundary particles in a coupled fashion (see Fig. 1 left flow chart).However, this coupling of boundary and fluid interaction causes an overestimation of SPH quantities for both interactions.More precisely, the density and, subsequently, the repulsive pressure force between the given fluid particle and the boundary depends on the number of fluid neighbors in its proximity.Thus, an individual fluid particle without any fluid neighbors will have less density and will receive less repulsive pressure force than a fluid particle with some fluid neighborhood.In practice, this yields a varying distance between fluid particles and the boundary, where fluid particles with many fluid particle neighbors have a larger distance to the boundary surface than a fluid particle with few or no fluid particle neighbors (see Fig. 2).Likewise, the calculation of fluid-fluid pressure or other forces, e.g., surface tension and artificial viscosity, within the boundary interaction radius are affected by the overestimation of the fluid particle's density due to boundary interaction, resulting in higher inter-fluid-particle forces, and, hence, in unintended fluid dynamics.initialize a f i by a i ▷ a i : non-press.acc.

13:
for all fluid particle i do update v i , update positions 31: end while was later improved by Band et al. 2018a.The concept of artificial density is additionally used to the regular density.Artificial density is dimensionless and represents the ratio between the current regular density ρ and the rest density ρ 0 , i.e., ρ a i = ρ/ρ 0 .Artificial density values ≥ 1, for example, result in an overcompression, yielding corresponding pressure forces.
First, we reformulate the artificial density calculation by separating it into three components: where i is the index of the current fluid particle, j is the index for all neighboring particles, f is the index for neighboring fluid particles, and b is the index for neighboring boundary particles.V i,j,f,b, denote the rest volumes of the particles (similar to Band et al. 2018a) and W ij is the kernel value between particles i and j (likewise for indices f, b).The first component in Eq. ( 1) is the "self-contribution" of the particle, the second component is a fluid neighbors' contribution, and the third component is a boundary neighbors' contribution.
For our method, we derive two densities.The fluid induced density for the fluid particle i is calculated by omitting the boundary neighbor contribution and the boundary induced density for the fluid particle i is calculated by omitting the fluid neighbors contribution component The main idea of this separation is that ρ f i represents the pure fluidfluid interaction which is not affected by neighboring boundary particles, whereas ρ b i represents the pure fluid-boundary interaction not affected by neighboring fluid particles.
Analogous to Band et al. 2018a, we use a constant boundary density ρ bb i = γ.However, as we use γ only for the boundary induced density without considering fluid particles, we use a slightly reduced boundary density of γ = 0.6 instead of γ = 0.7 as proposed by Band et al. 2018a.
Finally, the given iterative pressure solver is split into two stages, which are executed alternately, as shown in Alg. 1.The first stage comprises the calculation of the boundary pressure induced by fluidboundary interaction using ρ b i , while the second stage calculates the fluid pressure induced by fluid-fluid interaction using ρ f i .For both stages, the error thresholds, lim b and lim f , and the maximum number of iterations iter max are used to control convergence.These calculations are alternatively executed until both of the errors, err b and err f , are less than their thresholds or the number of global iterations iter g is more than iter max .Note, that the calculation of the current boundary pressure p b i requires the fluid induced acceleration a f i , and, correspondingly, the calculation of the current fluid pressure p f i requires the boundary induced acceleration a f i .The initial fluid acceleration is taken from the non-pressure forces.

Results and Evaluation
For the simulation, we use GPU-based opensource SPH framework openMaelstrom Winchenbach (2018).The radius of the fluid and boundary particles is set to 0.05 m which corresponds to approximately 0.1 m smoothing length.The constant γ for the density of the boundary particles is set to 0.6.For stabilization of fluid-fluid interaction, we use the artificial viscosity by Monaghan Monaghan (2005), with viscosity coefficient α = 0.05.We, furthermore, use the surface tension model by Akinci Akinci et al. (2013), with a tension coefficient of κ = 0.15.for the discretization and simulation, we use Cubic spline kernel Monaghan (2005) and Spiky kernel gradient Müller et al. (2004) for values and gradients, respectively.The gravity is set to (0, 0, −9.8) ⊤ m/s 2 .In the experiments with resting fluids in Secs.4.1, 4.2, and 4.5 we use regular rectangular samplings with varying spacing between boundary particles, while in experiments with fluid dynamics in Secs.4.3, 4.4, and 4.6, we randomly sample the boundary and optimize the boundary particle positions using the method from Bell et al. 2005 (with parameters: 100 impulse iterations, 0.01 impulse scale, 2 density).The rest density of the fluid is set to 1, 000 kg/m 3 .We use a fixed time step of 1 ms, which satisfies the CFL condition in all performed experiments.As the stopping limit for the fluid solvers, we use an average density error of 0.0001%.For evaluation, we implemented the pressure mirroring (PM) Akinci et al. (2012) and the pressure boundary (PB) Band et al. (2018a) methods, and we combined them with the IISPH solver in Band's modification Band et al. (2018a).To initialize the experiments, we relax the fluid using the relaxation stopping criterion by Akhunov et al. 2023.All evaluations have been performed on an AMD Ryzen Threadripper 3970X CPU with 64 GB of RAM and an NVIDIA GeForce RTX 4090 GPU with 24 GB of VRAM.
In all experiments, we denote the standard and decoupled versions of PB and PM as PBS, PBD, PMS and PMD, respectively.

Resting Fluid Sheet
In this experiment, we compare standard and decoupled versions of the PB and PM boundary handling methods, where the fluid sheet stays relaxed.The fluid sheet consists of 508 fluid particles lying on the regularly sampled boundary surface.The results of the experiment are presented in Fig. 3 and Tab. 1.We additionally vary the sampling distance between the boundary particles to investigate its influence on the particle height and stability.We use the term "sampling distance" as a ratio between the absolute distance between the neighboring boundary particles and the boundary particle diameter.We choose the sampling distances of 0.7, 0.6, and 0.5.For each distance, we measure the average height of the fluid over time and the average amplitude of velocity over time.
As can be seen from Fig. 3, the varying sampling distance changes both fluid height and fluid velocity for all methods.The larger the distance between boundary particles, the lower the height of the fluid.The decoupled versions of the methods, however, show a consistent height value for the PM and PB methods, while standard versions differ from each other significantly.Consequently, the stan- The cost for less noise and consistent fluidboundary distance is the number of iterations; see Tab. 1.For all sampling distances, the decoupled versions require more iterations than the standard, coupled versions, although the difference between the standard PM and the decoupled one is minimal.

Resting Fluid Bulk
This experiment is similar to the prior experiment, but instead of the sheet, a bulk fluid in a container is simulated.The fluid bulk consists of 4, 684 fluid par-  (2023).The 2D histograms indicate a significantly reduced level of noise for the decoupled versions of the boundary handling methods (see Fig. 5).There is a significant decrease in the residual pressure and the resid-  There is a significant difference in both, the pressure and the velocity distributions for PM, whose standard version is known to be very noisy Band et al. (2018a).The decoupled version of PM is less noisy compared to the standard one, and even less noisy than the decoupled version of PB.Fig. 4 visualizes the residual velocity amplitude for the fluid bulk experiment.It can be observed that the decoupling significantly reduces the residual velocity of the fluid particles close to the boundary and the free surface of the fluid.
As for the fluid sheet experiment in Sec.4.1, the number of iterations for this experiment is higher for the decoupled versions, as depicted in Tab. 2. For PM, the difference in iterations is small, similar to the prior experiment.

Fluid Sheet with Tangential Motion
In the next experiment, we compare the standard and the decoupled methods for a fluid sheet in tangential motion over the boundary surface.We sample the boundary surface randomly and then opti-mize the location of the boundary particles using the method by Bell et al. 2005.The fluid sheet is relaxed according to Akhunov et al. 2023, then accelerated in x-direction up to 1 m/s, and, finally, it slides freely without external acceleration besides gravity.In Fig. 6, we present the first 10 s of the sliding sheet.Similar to the experiment of the resting fluid sheet in Sec.4.1, the fluid height depends on the boundary handling method in the case of the standard version, while the height of decoupled versions is very similar.Over the time of 10 s, the average velocity decreases more quickly for the decoupled versions than for standard ones.This happens because of the reduced height of the fluid for the decoupled versions, i.e., the fluid lies closer to the boundary, and, hence, the boundary particle "features" cause more collisions with the fluid particles.These collisions dissipate the kinetic energy of the fluid particles.

Released Cylinder
Similar to the previous experiment in Sec.4.3, this experiment investigates the influence of the decoupling on a fluid's tangential motion.In this case, a fluid bulk is in a cylindrical container and after relaxation it is released, yielding more complex motion than in a single fluid sheet.The ground surface is sampled randomly and optimized with the method by Bell et al. 2005.The number of fluid particles is 103, 166.Due to the varying fluid-boundary distance, the height of the relaxed fluid depends on the bound- ary handling method and on the coupling/decoupling scheme applied; see Tab. 3, top row.Therefore, we enlarged the initial cylinder of 1 m radius to 3 m to achieve similar fluid heights, which leads to comparable potential energy in the system.In this case, the difference between heights is small, as it is shown in Tab. 3, bottom row.Fig. 7 depicts the top view of the scene 1 s after the fluid was released.It can be observed, that the advancing fluid front reaches slightly larger distances for the standard versions of the boundary handling methods than the distance reached by the decoupled version.This is due to the higher number of collisions

Water Jet Passing a Grid
In this experiment, a water jet is ejected on a grid with regular quadratic holes.The size of the holes is 2 × 2 particle diameter, and the rod width between the holes is also 2 particle diameter.The results are presented in Fig. 8.The results show that both standard coupled approaches, PMS and PBS, allow only a limited amount

Dragon Dam Break
The last experiment consists of a fluid volume of 745k fluid particles hitting a rigid dragon.The dragon is randomly sampled with 47k boundary particles and optimized with the method by Bell et al. 2005.The results of the experiment are presented in Fig. 9.
The main difference between the coupled and decoupled methods in Fig. 9 is that the coupled methods show much more splashy fluid behavior, which can be explained by unintended fluid motion near the boundary discussed in the "Resting Fluid Bulk" experiment in Sec.4.2.The overestimation of the fluid-boundary forces leads to the overestimation of fluid particles' velocity in the vicinity of the boundary and as a result brings an excess fluid motion, which in this case is a splashy fluid flow.
Some more subtle differences between the coupled and the decoupled method get visible around the fine details of the dragon.Some examples are shown in the insets in Fig. 9, showing close-ups of the dragon's tail (bottom left) and dragon's front claw (bottom right).These close-ups show that the decoupled density estimation yields a fluid flow that is significantly closer to the dragon in small cavities and that the flow bypasses the dragon's volume more smoothly.

Discussion and Limitations
Generalizability.While we demonstrated our decoupled boundary handling for particle-based bound-ary representations only, our decoupling method can be applied to any boundary handling approach which is based on density estimations.In the case of boundary integrals, for example, the iteration over the boundary particles for the density estimation can be replaced by the boundary integral calculation.Based on this density, the boundary pressure can be calculated using the boundary integral.In principle, the decoupling algorithm remains unchanged.
Performance.Tab.1 states that our method demands more iteration to converge than the regular methods.For the PB method, the difference is quite significant, however, for the PM method the difference is tolerable and only mildly affects the overall simulation time and performance.Depending on the application, our method may demand a trade-off between computation time and consistent fluid behavior.
Convergence.Although our method converges to the same value of the average density error, the residual fluid velocity in the vicinity of the boundary is much smaller, which makes the fluid behavior much more consistent between bulk regions and the regions near the boundary.
Fluid-boundary coupling.So far, our method supports only one-way coupling between fluid and boundary, and further investigations are required to combine our decoupling approach with methods like the Akinci et al. 2012or Gissler et al. 2019 to achieve two-way coupling.

Conclusion
In this paper, we propose a new method that decouples forces induced by fluid-fluid interactions and fluid-boundary interactions for SPH-based fluid simulations.The decoupling is realized by separately estimating the density and calculating pressure for fluidfluid and fluid-boundary interactions.This method prevents the overestimation of fluid density near the boundary and prevents too large distances between the fluid particles close to the boundary and between fluid and boundary particles.Preventing density overestimation makes the fluid flow more accurate and consistent with the bulk flow close to the boundary and at narrow and detailed boundary features.For example, the flow through small holes can be simulated accurately, while for standard methods small holes become almost impenetrable.Moreover, the decoupling reduces the amplitude of residual pressure and velocity, and it yields fewer oscillations at free surfaces.As a drawback, our method requires a larger number of iterations to converge and requires more investigation toward two-way coupling between fluid and boundary.

Figure 1 :
Figure 1: Depiction of the fluid and boundary neighbors of the selected (bold blue line) fluid particle.Arrows represent an interaction between the selected fluid particle and its neighboring fluid particles (blue) and boundary particles (yellow).Blue lines represent interaction with fluid neighbors, and yellow arrows represent interaction with boundary particles.The left flowchart is the classical approach in density and pressure calculation, and the right flowchart depicts our decoupled approach

Figure 2 :
Figure 2: Comparison of the fluid-boundary distance for a single fluid particle (left) and a fluid sheet with several fluid neighbors.The boundary particles and fluid particles are colored yellow and blue-purple, respectively.The top gray line shows the height of the fluid on the right.Screenshot from our SPH simulation while err f ≥ lim f ∨ iter ≤ iter max do 21:

Figure 3 :
Figure 3: Resting fluid sheet experiment.The average height values over time (top) and the average velocity amplitude (bottom) are shown.The experiment is done for different sampling distances.The first 10 s are shown after the relaxation criteria have been reached.

Figure 4 :
Figure 4: Resting fluid bulk experiment.Comparison of the pressure boundaries (PB * ) and pressure mirroring (PM * ) methods in their standard ( * S) with the proposed decoupled ( * D) implementations.The colormap encodes the velocity amplitude between 0 m/s (yellow) and 0.1 m/s (dark blue) 8 s after the fluid's relaxation is completed.

Figure 5 :
Figure 5: Resting fluid bulk experiment.The pressure and velocity histograms over time are presented.The measure of pressure and velocity is done for particles in bulk and for particles interacting with the boundary separately.

Figure 6 :
Figure 6: Fluid Sheet with tangential motion.The height values and the average velocity of the fluid particles during 8 s of sliding after the acceleration are shown.

Figure 7 :
Figure 7: Cylinder experiment.Top-down view of the distribution of the particles 1 s after cylinder release.

Figure 8 :
Figure 8: Water jet passing a grid experiment.The close-up of the jet going through the grid.

Table 1 :
No motion fluid sheet experiment.Number of solver's iterations

Table 2 :
No motion bulk experiment.

Table 3 :
Cylinder experiment.The height of the relaxed fluid bulk in a cylinder with a radius of 1 m and 3 m.