Monolithic Friction and Contact Handling for Rigid Bodies and Fluids Using SPH

We propose a novel monolithic pure SPH formulation to simulate fluids strongly coupled with rigid bodies. This includes fluid incompressibility, fluid–rigid interface handling and rigid–rigid contact handling with a viable implicit particle‐based dry friction formulation. The resulting global system is solved using a new accelerated solver implementation that outperforms existing fluid and coupled rigid–fluid simulation approaches. We compare results of our simulation method to analytical solutions, show performance evaluations of our solver and present a variety of new and challenging simulation scenarios.


Introduction
Rigid body dynamics, as well as fluid simulations are subject of an extensive amount of research motivated by great interest in applications in computer graphics, robotics and industrial prototyping [BET14,IOS*14,Bri15]. Even though most of the work focuses on exclusively simulating either the dynamics of rigid bodies or the motion of fluids, there exist more recent approaches combining fluids and rigid bodies into one simulation framework [AIA*12, MMCK14, KB17, HFG*18, GPB*19, BKWK19]. It is easy to see that these combined methods vastly expand the range of possible applications and simulation scenarios. Usually, in combined approaches, rigid contacts are resolved using a body representation and simulation method distinct from the fluid representation, fluidrigid interface handling and internal fluid pressure force computation [AIA*12, MMCK14, KB17, HFG*18]. Even in approaches that resolve rigid body contacts based on a particle representation, a concept well known from, e.g. Refs. [CS79,BYM05,TSIHK06,Ngu07,MMCK14,Coe17], the rigid-rigid contact handling typically still uses different methods than the fluid-rigid interface handling and the fluid solver. For instance, He et al. [HBH*18] and Peng et al. [PZWZ21] use smoothed particle hydrodynamics (SPH) to simulate fluids, yet employ a discrete element method to handle rigid contacts and fluid-rigid interaction. Similarly, Macklin et al. [MMCK14] implement a position-based SPH fluid, but on the other hand, a collision detection method similar to the molecular dynamics method [BYM05] together with shape matching constraints is used for rigid bodies. Due to impulsive collision responses, rigid contact handling strongly influences fluid-rigid constraints and consequently the internal fluid state. Equivalently, internal fluid forces might also have a great effect on neighbouring rigid bodies that needs to be considered when solving rigid contacts. By solving fluid and rigid body contact constraints sequentially, the mutual dependencies between the constraints are neglected, leading to unstable and inefficient coupled simulations. However, due to heterogeneous solving procedures for rigid bodies and fluids, mixed simulation methods have difficulties solving fluid and rigid body constraints in a combined manner. To overcome these issues, Gissler et al. [GPB*19] propose to simulate fluids, rigid bodies and their interactions exclusively with SPH and as such are able to build an intuitive and stable strong coupling between particle fluids and rigid bodies. However, even though realistic dry friction forces are crucial in order to authentically replicate rigid body behaviour, the dry friction force proposed by Gissler et al. [GPB*19] lacks basic features such as the ability to reproduce stiction and conservation of momentum. Additionally, despite the fact that it is common practice to compute contact and dry friction forces at once in order to maintain simulation stability [KSJP08,BET14], the explicit friction proposed by Gissler et al. [GPB*19] is applied separately from the contact handling procedure.

Figure 1:
A pure SPH simulation showcasing the capabilities of the proposed particle-based strong coupling between pressure forces applied in the fluid and at rigid bodies as well as dry friction forces between rigid bodies. Particles are used to represent the simulated rigid bodies, the shown fluid and boundary geometry. The tower, consisting of over 600 blocks, is brought to fall by two rigid ducks. When collapsing into a basin filled with fluid, our solver is able to robustly handle fluid incompressibility, fluid-rigid interaction, rigid-rigid collision forces and dry friction forces.

Contributions
Motivated by promising results of existing coupled simulation methods, our overall goal in this paper is to develop a monolithic simulation method that implements a strong coupling between fluids and rigid bodies which is as far-reaching as possible. For this, we build one global system that encodes all constraints for fluid volumes, fluid-rigid interface handling and rigid-rigid contacts. In contrast to existing coupled approaches such as the interleaved Jacobi solvers proposed by Gissler et al. [GPB*19], our unified system not only allows a more precise analysis of system properties (e.g. diagonal elements), it also enables us to use a non-linear conjugate gradient method, promising increased solver convergence speed. Since the capabilities of our rigid contact handling should match those of a designated rigid body simulator, we also develop a new implicit particle-based dry friction formulation for rigid-rigid contacts. This friction formulation addresses many shortcomings of the explicit friction force used by Gissler et al. [GPB*19], including the preservation of momentum and the ability to reproduce stiction. As the major improvement to Gissler et al. [GPB*19], our friction formulation is directly embedded in the global system and as such friction forces are included in the strong coupling between fluid pressure forces and rigid contact forces. To our best knowledge, there exist no other simulation method that implements such a far-reaching strong coupling between fluids and rigid bodies.
In summary, in this paper, we present • a novel particle-based implicit dry friction formulation based on the exact Coulomb friction model • a new monolithic, implicit SPH solver that unifies fluid pressure, fluid-rigid interface forces and rigid-rigid contact handling including dry friction forces into one system that is solved using a non-linear conjugate gradient method • results that validate our friction formulation and additionally demonstrate the capabilities and versatility of the proposed simulation method, as already hinted at in Figure 1.

Related Work
Computational methods to simulate rigid body contact dynamics have been studied extensively since quite a long time [MW88, Hah88, Bar89, Bar90, Bar91, Bar93a, Bar93b,Mir96,BB99,Ste00]. We refer to the report written by Bender et al. [BET14] for a great summary of earlier work and a general introduction to rigid body simulation. Previous work has shown that computing frictional forces is challenging due to their non-linear and nonsmooth relation to relative velocities and normal forces [ [YN06,XZB14]. In contrast to the methods described above, our friction computation is based on the exact Coulomb friction model including stiction, and we are able to directly embed the principle of maximum dissipation into our implementation without any further error-introducing simplifications. Newton-based approaches have shown to be able to implement the exact Coulomb friction model as well, by formulating the complementarity problem as non-smooth functions whose roots are found using a generalized version of Newton's method [BDCDA11, DBDB11, KTS*14, MEM*19]. Newton's method is quadratically convergent, which is especially useful when solving poorly conditioned problems, but in each iteration, a linear system needs to be solved [MEM*19, AE21]. The solver used in our implementation can also achieve superlinear convergence with only slightly increased computational costs per solver iteration compared to a standard Jacobi iteration [SHNE10]. Proximal operators, as already employed in Refs. [JM92,JAJ98], can be used as a tool to implement a friction force computation following the exact Coulomb friction law [Erl17]. Additionally, implementations based on proximal operators have shown to be easily extendable to employ other friction models such as anisotropic friction [EMAK19] and models considering frictional torque [LG03,Erl17]. Our proposed simulation approach utilizes proximal operators as the basis for a more advanced conjugate gradient method.
Even though implicit formulations have gotten more attention in the last few years, explicit pressure solvers are still widely used in the context of SPH fluids [BT07, AIA*12, IOS*14, RHEW17]. For less complex scenarios (e.g. low fluid depth, softer requirements regarding incompressibility), explicit formulations benefit from lower computational costs per simulation step while the increased simulation stability through implicit formulations cannot outweigh their higher computational burden [KBST19,KBST22]. Similarly, explicit penalty methods for rigid body simulations are easy to implement and require little computational cost per step [MW88,BET14]. For our simulation method however, we still chose to employ an implicit formulation. This not only allows us to handle rigid contact constraints similar to fluid incompressibility, we are also able to simultaneously solve for pressure in the fluid, pressure forces at the fluid-rigid interface and between rigid bodies as well as friction. Computing friction and contact forces at once is known to be beneficial for simulation robustness and performance, since frictional forces can have a significant effect on pressure constraints at contacts and vice versa [KSJP08,BET14]. In Section 5.1, we demonstrate that explicit friction formulations, such as the one proposed by Gissler [BT07] or PCISPH [SP09] and deriving pressure forces that are then applied to fluid and rigid particles. Friction is only mentioned to be computed at the fluid-rigid interface employing a laminar artificial viscosity force. Rigid-rigid contacts, including dry frictional forces, are handled separately using an external simulation software such as Bullet [Cou]. In contrast, our simulation method unifies fluid pressure forces with rigid contact handling and friction, resulting in more stable fluid-rigid interfaces as already demonstrated in Gissler et al. [GPB*19]. Macklin et al. [MMCK14] introduce a position-based particle simulation framework that combines many object types, including fluids and rigid bodies. Particles representing rigid bodies are treated as if they were unconnected during a solver iteration and shape matching constrains are applied afterwards to ensure rigidity. It is mentioned that the cost of shape matching grows quickly with the number of particles and impulse propagation is rather slow [MMC*20]. They present a pairwise particle friction model for granular materials which is also applied to rigid body particles. Frictional forces are computed and applied iteratively during the constraint solving procedure of a simulation step. The authors state that the simulation method aims at real-time performance while making compromises in realism, and as such there is no validation shown for the correctness of frictional forces depending on the parameterizable coefficient of friction. Instead, the authors report that frictional forces strongly depend on the iteration count of the constraint solver [MMCK14]. The development of extended position-based dynamics (XPBD) [MMC16] mitigates the dependency of constraint forces on the iteration count in position-based simulation approaches but does not explicitly discuss frictional forces. Subsequent publications based on XPBD do not demonstrate the physical correctness of the computed frictional forces either [MMC*20]. Our simulation method uses the coefficient of friction μ as a physically meaningful input parameter and is able to replicate the expected friction behaviour. In their coupled fluid-rigid simulation approach, Koschier and Bender [KB17] externalize the rigid-rigid contact handling, but still apply an explicit Coulomb friction force on fluid particles next to a rigid boundary that uses pressure values to estimate normal force magnitudes. While our implicit friction implementation will use a similar idea to estimate normal force magnitudes between rigid bodies, we also determine contact normal directions based on SPH pressure forces. This way, in contrast to Koschier and Bender [KB17], no boundary geometry description other than the particle representation is required. Hu et al. [HFG*18] present an MPM method implementing a two-way coupling between rigid bodies and objects such as fluids, elastic and elastoplastic materials. However, for rigid-rigid contacts they also use an external rigid body dynamics software. Accordingly, frictional forces between rigid bodies are not discussed. In their partitioned approach, Akbay et al. [ANZS18] couple blackbox solvers for fluids, other deformable objects as well as rigid bodies through a small reduced-order system. As an example, this way they are able to couple the Lagrangian fluid solver DFSPH [BK15] with the position-based rigid body solver presented by Deul et al. [DCB14]. Partitioned simulation methods have the advantage that existing specialized solvers can be reused for a combined simulation. However, the authors mention that in general, compared to partitioned approaches, monolithic approaches that unify all simulated quantities in one system are expected to be more efficient in computing a strongly coupled solution [ANZS18].
More recently, Gissler et al. [GPB*19] presented a partitioned simulation method to couple particle fluids with particle-based rigid bodies. Both, the fluid solver and the rigid body solver are purely based on SPH, such that interleaving both solvers is particularly straightforward. The rigid body solver detects potential collisions between bodies by considering density deviations at particles and prevents interpenetrations using pressure forces. To model friction between rigid bodies, they propose to use an explicit Coulomb friction implementation. As we will demonstrate in Section 5.1, this friction implementation cannot reproduce static friction and is not guaranteed to result in the physically correct amount of friction force. Also, frictional forces are not guaranteed to preserve momentum. Our rigid body solver uses density deviations to detect contacts and pressure to resolve collisions as well, however, in contrast to Gissler et al. [GPB*19], we present a truly monolithic solver that strongly couples rigid body contact handling with a particle fluid simulation. Additionally, the pressure system is extended by a novel momentum-conserving implicit Coulomb friction force formulation that correctly reproduces static friction and stick-slip transitions. This way, we need to solve only one global system to simultaneously handle fluid incompressibility constraints, fluid-rigid interface forces, rigid-rigid contacts as well as rigid-rigid friction forces. The unified system makes it easier to use more advanced solving methods, such as the non-smooth non-linear conjugate gradient method (NNCG) proposed by Silcowitz et al. [SHNE10]. Compared to the relaxed Jacobi method employed by Gissler et al. [GPB*19], we are able to show performance gains above a factor of 10.

Method
A complete rigid body contact handling procedure includes resolving potential collisions between simulated bodies and computing corresponding frictional effects. In this section, we describe the main concepts of our monolithic simulation method that is purely based on SPH discretizations. Section 3.1 shows how SPH can be used to detect contacts between rigid bodies as well as contacts at the interface between fluid and rigid bodies. The definition of contact is then used in Section 3.2 to construct the pressure LCP which forms the main optimization problem that needs to be solved in order to handle rigid contacts and fluid incompressibility. Similarly, in Section 3.3, we introduce the optimization problem for frictional forces following the exact Coulomb friction model. We will see that both optimization problems form a system since they mutually depend on each other. Afterwards, in Section 4, we describe how this system can be solved in an implementation.

Contact detection
Typical simulation methods, that are handling contacts between bodies represented as triangle meshes, perform a collision detection step at the beginning of each simulation iteration in order to identify all contact points in the scene. This is a non-trivial task and the quality of the contacts generated by the collision detection can severely impact performance, robustness and correctness of the simulation result [BET14, Erl18, AE21, WFS*21]. To make matters worse, good quality of contact points is elusive and neither rigorously understood nor defined [Erl18]. Additionally, comparing simulators becomes difficult as often the underlying contact generation method is not sufficiently specified [Erl18]. In our particle-based simulation framework on the other hand, inspired by previous work (e.g. Refs.[AIA*12, GPB*19]), the surface of rigid bodies is sampled with SPH particles such that collision detection can be boiled down to a volume estimation of particles using a local particleneighbourhood. Neighbouring particles can be reliably found using well-studied neighbourhood search methods [IABT11, IOS*14, WSG*18, BGT19]. Using SPH, we compute the volume V r of a particle r that lies on the surface of a rigid body R with where r r are neighbouring rigid particles belonging to the same rigid body R as r, r k represent neighbouring rigid particles belonging to other rigid bodies and r b are neighbouring particles that are part of a kinematic boundary. W is the SPH smoothing kernel function. We use the notation W rrr shorthand for W (| x r − x rr |) with particle positions x. Additionally, if no time index is given we refer to current positions x(t , the rest volume V 0 r of particle r is approximated at the start of the simulation by considering surrounding particles r r sampling the same rigid body R: We now define that a particle r belonging to a rigid body is in a state of collision if In Equation (2), γ is a correction coefficient that accounts for the fact that rigid bodies are only sampled at the surface, whereas an SPH approximation of the volume assumes a complete particleneighbourhood. With γ , we can control how many contributions from neighbouring particles are necessary for a particle r to be considered in a state of collision [GPB*19] as illustrated in Figure 2. In our simulations, γ is set equal to 0.7 as we found this value to be high enough such that contacts between small and thin structures are still detected.
Traditionally, collision detection algorithms search for vertexface and edge-edge collisions with the aim of returning a set of collisions that individually describe the contact between exactly two bodies [BFA02, KSJP08, OTSG09, BET14, AE21]. As a major difference to most existing approaches, we do not explicitly define contact points between rigid bodies since the volume of a particle r can be influenced by particles r k belonging to multiple rigid bodies nearby, and a contact between rigid bodies can cause multiple particles on both sides of the contact to be in a state of collision. Consequently, our contact handling and friction force computation will also be based on the particle description of rigid bodies instead of generated collision points.

Figure 2:
An illustration of the effect of γ on the contact detection between two rigid bodies. With γ = 0.3, interpenetrations between bodies may occur since contributions from neighbouring particles might not be enough to detect a collision. This is especially problematic for thin geometry features such as the yellow rod displayed. Setting γ = 1.0 causes large gaps between rigid bodies and finer features such as the notch in the grey body is ignored due to the collision between the bodies being detected too early. We achieve the desired result by setting γ around 0.7.
We want to point out that most recent SPH-based pressure solvers use a similar constraint on the volumes of fluid particles f to define incompressibility. Here, the solvers preserve a compression free fluid state where no fluid particle f has throughout the simulation to prevent any volume loss in the fluid [SP09, ICS*14, BK15]. Fluids are sampled volumetrically, so the rest volume V 0 f of fluid particles f is set to where h is the particle spacing. The current volume V f can be approximated using with neighbouring fluid particles f f , neighbouring particles f k that sample rigid bodies and neighbouring particles f b belonging to a kinematic boundary. In Section 3.2, we will see that the similarity between incompressibility constraints allows us to comfortably combine the volume constraints of rigid and fluid particles into one global system.

Discussion
Using particles to detect and handle rigid body collisions entails interesting differences compared to traditional mesh-based schemes. Precision of contact handling no longer purely depends on mesh geometry, but instead the sampling density of particles governs the resolution of contact handling. Existing work such as Allard et al.
[AFC*10] indicate that it can be advantageous to abstract contacts away from traditional vertex-triangle and edge-edge collision pairs. In the specific case of Allard et al. [AFC*10], one intersection volume is constructed per contact which is then sub-divided by a regular grid into multiple sub-volumes to resolve contacts more accurately. By adjusting the grid resolution and particle sampling density, respectively, Allard et al. [AFC*10] and our approach both can tune contact resolution independently of the complexity of the underlying mesh. As a consequence demonstrated in Section 5.4, highly detailed meshes can be simulated using a more appropriate precision. On the other hand, even though at first glance, it seems excessive to sample simple meshes with a high number of particles, the increased contact resolution allows us to simulate pressure and friction distributions over contacting surface patches since potentially all particles within the contact area are detected to be in a state of collision. This contrasts with most rigid body simulation methods which only consider a few distinct collision points per contact area [BET14,Erl18]. We refer to Allard et al. [AFC*10] for an interesting demonstration on how the number of contact points per colliding bodies is indeed relevant for replicating realistic body motion.
As a disadvantage of the employed contact detection using particle volumes, it shall be mentioned that if geometrically highly accurate collision points based on mesh descriptions of bodies are required, the number of particles necessary to achieve such a high resolution might significantly reduce simulation performance.
Parameter γ γ γ : Next to the particle resolution, parameter γ ∈ (0, 1) influences contact detection for rigid bodies by scaling the rest volumes of rigid particles. The effect of γ is illustrated in Figure 2. As already mentioned, decreasing γ towards zero increases the required contributions from neighbouring particles before a rigid particle is considered to be in a state of collision. Thus, we need to make sure that γ is set sufficiently large such that no thin body features, sampled with few particles, can penetrate through another body surface without causing a collision to be detected. On the flip side, setting γ too close to one causes particles to be already considered in collision even if rigid bodies are still quite far away from each other. This problem could be mitigated by offsetting particles along the surface normal towards the inside of the body. However, similar to Gissler et al. [GPB*19], as a more general solution requiring less hand tuning and that works for open meshes, we found that using γ = 0.7 results in a robust and accurate contact detection with no need to offset particles. A more theoretical motivation behind γ = 0.7 is given by Gissler

Pressure system
In this section, we define the pressure LCP that constrains all particle volumes. By solving the pressure LCP as shown in Section 4, we obtain pressure forces that preserve the volume of all particles in the simulation and as such we ensure that fluid incompressibility is met and no rigid bodies are intersecting. After defining the pressure LCP, in this section, we also derive the relations between pressure and particle volume since they are required in a solver implementation.
To derive the pressure LCP, we start by defining a volume error V err i at each particle i with and reconsider the incompressibility constraint from Equations (3) and (4) at the next timestep t + t for some fluid or rigid particle i: Note that only negative volume error corresponds to a compression of a particle, while positive volume error indicates that a particle is in a stress-free state with larger volume compared to its rest volume.
The goal of the pressure solver is to compute pressure values p i such that Equation (8b) holds true, meaning that the predicted volume error V err i at the next timestep t + t must be greater or equal to zero for all rigid and fluid particles i. Since we consider volume errors at particles at the next timestep V err i (t + t ), which already depend on current pressure, we obtain an implicit pressure formulation. Similar to Gissler et al. [GPB*19], we solve for unknown pressure values p i at rigid and fluid particles i and derive pressure forces F P i based on those pressure values. To prevent unphysical behaviour, we further restrict p i to be non-negative. This is a known method used in SPH pressure solvers [ICS*14, BK15] to prevent attraction forces between particles due to negative pressure and can be compared to the non-negativity constraint of normal forces acting at rigid body contacts [Erl07, Erl13, BET14, MEM*19, AE21]. Additionally, we only allow p i to be non-zero if the predicted volume error at the next timestep V err i (t + t ) equals zero. Together, the constraints on V err i (t + t ) and p i form an LCP [ANE17] which can be written down using the more compact notation Note that on a high level, this pressure LCP closely resembles LCPs known from traditional mesh-based contact handling methods [BET14, Erl13, MEM*19]. Typically, solving the LCPs means finding scalar values for normal forces, or in our case pressure values, that cause some collision term, often a relative velocity [BET14] but in our case V err i (t + t ), to become non-negative. Typical velocity-level contact handling approaches require constraint stabilization methods to prevent bodies from drifting into each other over time [TBV12,BET14,Erl17,AE21]. These methods often introduce ad-hoc forces into the simulation and are known for causing stability issues [BET14, Erl07, MEM*19]. In contrast, our simulation method explicitly computes volumes in each simulation step based on the distances to neighbouring particles and thus the necessity for any constraint stabilization that prevents intersections over time is eliminated. Due to the similarity of the LCPs, in principle, the same solving mechanisms employed in traditional mesh-based simulation methods can be used to solve Equation (9).
To be able to find a solution to Equation (9), we must know how V err i (t + t ) depends on all unknown pressure values p. Thus, in the following, we derive the relation between V err i (t + t ) and p for rigid particles r and fluid particles f .

Rigid particles:
We begin with the definition of the volume error from Equation (7) and also consider Equation (1) describing the computation of V r . Distances between particles r and r r belonging to the same rigid body never change, so contributions from the first sum are constant in time. In our simulation, we employ an Euler-Cromer time integration scheme. Together with a first-order Taylor expansion of W (t + t ), we can write Equation (10) approximates the dependency of the predicted volume error V err r (t + t ) on the velocity field v at the next timestep t + t. Velocities of kinematic boundary particles r b are pre-defined and not influenced by pressure forces. What remains is finding the relation between the velocity v r (t + t ) of a rigid particle r, pressure forces F P and pressure values p. We predict Here, v R and ω R are translational and angular velocities of rigid body R, R also represents the set of all rigid particlesr sampling R, M R is the mass of body R, x R is the centre of mass of R and I −1 R denotes the inverted inertia tensor of R. The velocities v * R and ω * R represent intermediate translational and angular velocity body R has right before the pressure solve. Thus, v * R and ω * R include all explicitly computed forces and accelerations which, for a rigid body, typically are gravity g and the gyroscopic force [Ben07,BET14]: Now, the last unknown in Equation (11) are pressure forces F P and frictional effects F F and τ F . Section 3.3 describes in detail how frictional forces are computed, however, their embedding in the pressure system already indicates the employed strong coupling between frictional and pressure forces. To compute F P , we use a symmetric SPH pressure gradient estimation that uses a volume formulation (as shown before by e.g. Band Our method is flexible about the computation of pressure values p r b at kinematic boundary particles r b , in the sense that any boundary handling method that defines a pressure value at neighbouring boundary particles r b (e.g. Refs. [AIA*12, BGI*18, BGPT18]) could be used to define p r b . In the following, we will assume that pressure is mirrored with p r b := p r . Also note that pressure values of neighbouring fluid particles p r f are considered in the pressure force estimation, enabling our strong coupling between rigid bodies and fluids. In summary, Equations (10) to (13) describe the relation between the volume error V err r of a rigid particle r and the vector of all unknown pressure values p.
Fluid particles: In our simulation method, pressure at fluid particles f and pressure at rigid particles r are computed simultaneously. Since both have the same constraint on the volume to enforce incompressibility, we can derive the dependency of V err f (t + t ) on p for fluid particles f in a very similar manner. Expanding Equation (8b) yields The computation of v f k (t + t ) for rigid particles f k is already known from Equation (11), the velocity where m f = V 0 f ρ 0 f is the mass of fluid particle f with fluid rest density ρ 0 . Similar to above, we use pressure mirroring at the boundary with p f b := p f . Due to their symmetry, Equations (13) and (15) guarantee that pressure forces conserve momentum exactly. As a small but important detail, we want to point out that pressure forces between fluid particle f and neighbouring rigid particle f k are solely defined by fluid pressure p f . We do not consider p f k to prevent highpressure accelerations of relatively lightweight fluid particles f next to a rigid-rigid contact that typically causes much higher but punctual pressure at contacting rigid particles f k . Again, v * f is the velocity of particle f including all forces computed prior to the pressure solve, so we have where F E f are additionally computed forces such as viscosity (e.g. [Mon92,WKBB18]) and surface tension (e.g. [AAT13]). Using Equations (10) to (16), the pressure solver presented in Section 4 aims at computing pressure values that fulfil Equation (9) at all fluid and rigid particles, and as such all at once prevents rigid bodies from penetrating into each other, handles interface forces between fluids and rigid bodies and guarantees fluid incompressibility.

Frictional forces
In the previous section, we defined the optimization problem that can to be solved for pressure values p that prevent compressions at all particles. Here, we do the same for friction forces. We first define the underlying optimization problem and then provide details that will be required by the solving procedure described in Section 4. We chose to model frictional forces according to the exact Coulomb friction constraint. More precisely, we are searching for friction multipliers λ r that are bound by the friction cone F r : where n r is the contact normal direction, μ r is the coefficient of friction and F N r represents some measurement of the normal force acting at particle r. Note that we distinguish between frictional multipliers λ and frictional forces F F or frictional torque τ F as friction multipliers λ r only consider frictional effects at particle r due to normal forces induced by pressure at particle r. Those frictional effects encoded in λ r need to be mirrored onto neighbouring rigid particles to result in physically correct frictional forces F F r and frictional torque τ F r . We will describe this in greater detail later in this section. In addition to the Coulomb constraint, frictional multipliers should follow the principle of maximum dissipation stating that from all possible λ inside the friction cone F r , the friction multiplier λ r should fulfil [Erl17] λ r = arg min λ∈Fr λ · v tang with tangential velocities v tang r . By considering tangential velocities v tang r (t + t ) at the next timestep t + t, we obtain an implicit friction force formulation.
Particle framework: Computing frictional forces in a particle framework imposes further challenges concerning the contact geometry. While mesh-based collision detection usually returns one contact point and a contact normal direction [AE21], for a particle that is detected to be in a state of collision, these properties need to be estimated from neighbouring particles. Gissler et al. [GPB*19] propose to use an SPH summation to estimate the contact normal direction n r at a particle r: Here, n solely depends on neighbouring particle positions and the normal force F N r needs to be estimated separately using collision impulses. Since our simulation method solves pressure and frictional forces simultaneously, we propose to directly use intermediate pressure values p to estimate F N r and n r . For this, we consider the contribution of p r to its pressure force F P r for each particle r as described in Equation (13): This does not only eliminate the need for a separate contact normal estimation process including all associated difficulties, but we can also directly relate pressure forces with the normal force used in the friction computation. In Section 5.2, we show that this new approach does not suffer from problems associated with wrongly estimated contact normals and thus trivially handles all degenerated test cases proposed by Erleben [Erl18]. We note that Koschier and Bender [KB17] as well as Winchenbach et al.
[WAK20] also use the pressure force magnitude as an estimate for the normal force magnitude when computing frictional forces at the fluid-rigid interface. However, they do not compute the contact normal directions based on normal forces but instead derive them from the boundary geometry.
Tangential velocities: Equation (18) requires the computation of tangential velocities v tang r (t + t ). They can be estimated with a standard SPH approximation to compute a relative velocity v rel r (t + t ) at particle r and projecting v rel r (t + t ) onto the tangential contact plane afterwards using n r : v rel with the identity matrix . From Equation (14), we already know how to compute v r (t + t ) and v r k (t + t ). v r b (t + t ) is independent of pressure and friction forces, so in this context, we set v r b (t + t ) = v r b (t ).
Friction force mirroring: Friction multipliers λ r need to be mirrored onto neighbouring rigid bodies to result in physically consistent friction forces. Equivalent to pressure forces, we distribute the friction multipliers considering the weighting by the kernel gradient. The actual friction forces F F r and torques τ F r applied to rigid particles r are then given by It is easy to see that the friction force F F r includes its own friction multiplier λ r as well as mirrored contributions from neighbouring rigid particles r k . The friction force mirroring guarantees exact conservation of linear and angular momentum. Even though momentum conservation is no concern when interacting with kinematic boundaries, neighbouring particles r b belonging to such a boundary are treated consistently. Equations (22a) and (22b) describe linear relationships between F F r , τ F r and the vector of all unknown friction multipliers λ, and as such can be smoothly embedded into our system relating v tang r (t + t ) with λ. The modified volume V r is easily computed on-the-fly as described in Equation (23). In Section 5.1, we demonstrate the correctness and accuracy of the computed friction forces.

Summary:
Since we employ an implicit force computation model, we need to be able to describe the relation between the vector of all unknown friction multipliers λ and v tang r (t + t ) for all particles r. In the following, we give a short summary on how v tang r (t + t ) is computed. Given the frictional multipliers λ specified by Equations (17) and (18), we start by mirroring λ onto neighbouring bodies using Equations (22a) and (22b), to obtain frictional forces F F and torques τ F for all rigid particles. Frictional forces and torques, together with pressure forces F P , can be plugged into Equation (11) to obtain an estimate of v r (t + t ) for all r. The newly estimated velocities are inserted into Equation (21) to compute v tang r (t + t ).
We point out that the estimate for V err i (t + t ) required in the pressure system as well as the estimate for v tang r (t + t ) required in the friction system both simultaneously consider pressure forces acting in the fluid, at fluid-rigid interfaces and between rigid bodies together with frictional forces. This way, we achieve a strong coupling between contact forces and friction through one global monolithic system, which is crucial for simulation stability when considering rigid bodies [KSJP08]. While this is common practice in pure rigid body simulation frameworks [BET14, Erl17, PAK*19, MEM*19, MEM*20, FLS*21], we are not aware of an existing simulation method that strongly couples internal fluid pressure forces, fluid-rigid interface forces, rigid-rigid contact forces with dry friction forces. In Section 4, we provide an illustration of an efficient solver implementation that is able to compute p and λ satisfying the constraints given in Equations (9), (17) and (18).

Implementation
The optimization problems presented in Equations (9), (17) and (18) form the foundation of the pressure and friction solving procedure. In this section, we now build a solver implementation that is able to efficiently solve the previously presented optimization problems for pressure p and friction λ. For this purpose, both, the pressure LCP (Equation (9)) and the friction computation problem (Equations (17) and (18)) are translated into an equivalent fixed point formulation [JAJ98, SNT11, Erl17]: where P i is the allowed solution space of p i and F r the set of all friction forces inside the friction cone: The proximal operator prox S projects its input onto the nearest point in some given set S: We refer to Schindler et al. [SNT11] for a thorough derivation of Equation (24) from the pressure LCP and the friction constraints. In Equation (24), the scalar parameters α P and α F have no effect on the analytical solution of a fixed point as long as α P , α F > 0. In an implementation however, they greatly influence convergence speed and stability of an iterative solver. By iteratively updating p i and λ r using Equation (24), one generates a sequence of iterates that converges locally if the spectral radius of the Jacobian of the right-hand side of Equation (24) is smaller than one [Erl17]. This requirement can be met by setting α P and α F sufficiently small [FGNU06,PB14]. While there exist approaches that employ backtracking strategies to find safe values for α P and α F [Erl17], we found that using values that correspond to a modified Jacobi update similar to Tonge et al. [TBV12] and Gissler et al. [GPB*19] also leads to a robust and performant simulation.

Jacobi scheme
We can now build an iterative Jacobi scheme to solve the fixed point problem in Equation (24)  24a) corresponds to a relaxed Jacobi update step with relaxation coefficient υ where for rigid particles, the number of contacts n R of body R is additionally taken into account: The terms ∂/∂ p i V err i (t + t ) are an inherent part of the Jacobi update. They are used as an estimate of how much V err i (t + t ) changes depending on the unknowns p i . Contrary to the suggestion of Erleben [Erl17], we do not use a Gauss-Seidel variant to solve Equations (24a) and (24b). While the increased stability might outweigh the loss of parallel computation power for meshbased rigid body simulators, in our particle-based contact handling method including fluids, there are easily hundreds of thousands of simultaneously active contacts at each simulation step.
Handling this large number of contacts requires a parallelized solving procedure. It has been shown that Jacobi methods require additional stabilizing when treating rigid contacts [TBV12,GPB*19], which is why we divide the relaxation coefficient by n R . The number of contacts n R is updated in each Jacobi iteration using where V err r (t + t ) is the current prediction for the volume error at particler at the next timestep. The interlinked pressure computation method by Gissler et al. [GPB*19] uses a generic pressure solver to compute pressure in fluids. However, generic fluid solvers are not able to capture the effect velocity changes of rigid particles f k due to pressure p f have onto the volume error V err f . Thus, there are missing contributions in the computation of ∂/∂ p f V err f (t + t ) that can impact the convergence behaviour of the Jacobi solver. In contrast, our method unifies the pressure computation for rigid and fluid particles into one system, and thus the dependency of v f k (t + t ) on p f can be correctly included in the computation of ∂/∂ p f V err f (t + t ). The derivation and an efficient computation implementation of the elements ∂/∂ p i V err i (t + t ) are shown in Appendix A.
Friction update: When applying a fixed point iteration for frictional forces as shown in Equation (24b), we found it to be important that α F truly is a positive scalar value. Blocked Jacobi schemes that use a matrix for α F [CPS09] often cause α F r v tang r (t + t ) to point into a different direction than v tang r (t + t ). This has the effect that frictional multipliers λ r fulfilling the principle of maximum dissipation are no longer a solution to the fixed point problem described in Equation (24b). Erleben [Erl17] evaluated the convergence behaviour of blocked solver strategies in a very similar context and the results seem to support our claim. Similar to the blocked approach, considering α F independently for each coordinate direction of λ r introduces the same problem, which is why in our implementation, we choose where tr is the trace operator. Again, an efficient implementation to compute ∂/∂ λ r v tang Further, we slightly modify the Jacobi update for frictional multipliers λ r given in Equation (24b). Based on the observation that for two rigid particles r 1 and r 2 with v tang λ r 1 and λ r 2 should converge to the same vector. However, by looking at Equation (24b), we notice that since v tang r 2 (t + t ) is much larger than v tang r 1 (t + t ), λ r 2 will converge much faster compared to λ r 1 . This causes severe stability issues in the solving process for large relaxation coefficients υ or alternatively impractically slow convergence for smaller υ. To stabilize the solving process, we introduce a friction target λ * r into the solving process whose magnitude is clamped at μ r | F N r | and reformulate the fixed point problem from Equation (24b) to Again, this does not change the analytical solution of the fixed point problem but merely improves convergence behaviour.

Non-smooth non-linear conjugate gradient
Relaxation methods such as the Jacobi method are known for their simplicity, robustness and flexibility, but also come with slow convergence speed, especially considering poorly conditioned problems [CPS09,SNE09,PAE10,Erl13]. To improve the performance of the proximal operator method, Erleben [Erl17] suggests studying a combination of proximal operators, embedded into a generalized conjugate gradient method. We follow this suggestion by extending our Jacobi iteration into a NNCG method as described by Silcowitz-Hansen et al. [SHNE10]. Here, differences between iterates of the projected Jacobi solver are seen as residuals r where PJ returns the iterates after applying one projected Jacobi iteration as described in Equations (24a) and (31). The residuals simultaneously equal the negative gradient of some function f (r k ) := 1 2 |r k | 2 . Finding p and λ such that f becomes zero is equivalent to solving the fixed point problems given in Equations (24a) and (31). To find a local minimum of f , the Fletcher-Reeves non-linear conjugate gradient method can be employed [NW06,SHNE10,AE21] as it only requires information about the gradient ∇ f = −r. The implementation of the NNCG method is illustrated in detail in Section 4.3. Performance comparisons between the standard Jacobi method and NNCG are shown in Section 5.3.

Algorithm
We use this section to give an overall view of our simulation method and discuss some remaining implementation details. Algorithm 1 shows a summary of our particle-based fluid-rigid simulation loop, which closely resembles the typical structure of a simulation based on SPH [IOS*14]. As we can see, a simulation step starts by searching and storing particle neighbourhoods as these are required multiple times later on in various SPH summations. Similarly, the current volume V i (t ) is pre-computed for all rigid and fluid particles i using Equations (1) and (6), as its value is repeatedly used. Using SPH, we are able to explicitly compute forces such as viscosity and sur-Algorithm 1. A fluid and rigid body simulation loop based on SPH particles. The pressure and friction solver embedded in line 1 is presented in detail in Algorithm 2.
face tension and apply them using Equation (16) (12). The intermediate velocities are used in the pressure and friction solver which is embedded into the simulation loop in line 1. For now, we treat the solver as a black box and end the simulation step by integrating the resulting pressure and frictional forces onto fluid particles f , rigid bodies R and rigid particles r. A precise description of how velocities and positions are integrated in time is moved to Appendix B as the expressions become quite lengthy without adding much to the understanding of the simulation method.
The pressure and friction solver shown in Algorithm 2 is the core of our simulation method. Based on predicted particle velocities v i (t + t ) we implicitly compute pressure values that prevent collisions and matching frictional forces. The implemented NNCG solver bases on Jacobi updates which are described in Algorithm 3. Before each Jacobi step, we need to evaluate predicted volume errors V err i (t + t ) and tangential velocities v tang r (t + t ) based on current pressure vector p and friction vector λ. While this procedure is computationally expensive, we are able to implement it in an efficient and fully parallelized manner as demonstrated in Algorithm 4. To make Algorithm 2 easy to reimplement, a detailed explanation on how to compute the diagonal elements is given in Appendix A.
Typical SPH pressure solvers use the average predicted volume error of particles as a measurement for convergence [ Algorithm 2. A NNCG implementation to solve for unknown p and λ using the projected Jacobi step as an estimate for the function gradient. Afterwards, we can compute the desired pressure forces F P i and frictional effects F F r and τ F r from p and λ.
Algorithm 3. The projected relaxed Jacobi step to update pressure p and friction multipliers λ. Note that in line 3, we use the improved Jacobi step from Equation (31) to update λ r .

Results
In this section, we demonstrate the capabilities of our monolithic solver. We begin with a validation of our friction computation method by comparing simulation results to analytical solutions in Section 5.1 and perform some basic robustness tests in Section 5.2.
To motivate the usage of the NNCG solver, we compare its convergence speed to that of a standard projected relaxed Jacobi solver in Section 5.3. More challenging scenarios are showcased in Section 5.4, including difficult frictional settings, large mass ratios, complex geometries, interaction with fluids and large numbers of simulated rigid bodies. Please refer to the accompanying video for a more detailed insight into the simulation scenarios.
Algorithm 4. The procedure to compute predicted volume errors V err (t + t) and tangential velocities v tang (t + t) based on current pressure values p k and friction multipliers λ k . Since predicted velocities of rigid particles v r (t + t) can be computed on-the-fly based on v R (t + t) and ω R (t + t) using Equation (11a), we only have to iterate twice over all particles.

Validation
We use a rigid body placed onto an inclined plane as a basic test for the precision of our friction formulation. Given the slope of the plane θ, it is known that the rigid body R should slide down the ramp if μ R < tan θ with coefficient of friction μ R . In our test scenario, the plane has a slope of θ = 30 • , we start with μ R = 1 and reduce μ R with constant speed over time. Figure 3 visualizes the scenario and gives a first impression of the test results. As we can see, the friction model proposed by Gissler et al. [GPB*19] is not able to replicate stiction as the rigid body slides down the slope even for μ R = 1. On the other hand, our friction formulation correctly causes the body to stick to the surface as long as μ R > tan 30 • ≈ 0.577. Figure 4 gives a more detailed insight into the relation between coefficient of friction and velocity of the rigid body simulated using our friction formulation.

Solver robustness
To hint at the robustness of our particle-based contact handling, we simulate the degenerated test cases proposed by Erleben in Ref. [Erl18]. These special-cases of contacts between bodies are challenging for traditional mesh-based contact handling schemes [FLS*21]. As shown in Figure 5, our simulation method does not show any difficulties when simulating these constellations. However, we still like to mention that contact normals computed in our simulation only approximately match the ones listed by Erleben in Ref. [Erl18] due to the particle sampling of surfaces. The small errors in the contact normal computation have negligible effect on the overall simulation results which is why we do not consider them any further, but instead refer to work published by Koschier and Bender [KB17] as well as Bender et al. [BKWK19] that addresses a very similar issue.

Figure 4:
The relation between the coefficient of friction μ R of rigid body R and the magnitude of its translational velocity v R using our friction model. We can see that R starts accelerating as soon as μ R < 0.577 which is in good agreement with the analytical solution. Note that the velocity increases in a quadratic manner as μ R decreases.

Solver performance
Large mass ratios between interacting bodies are traditionally hard to handle. As shown in Figure 6, our method is able to stably simulate a rigid cube resting on top of a second cube, whereby the upper cube is 1000 times heavier. We use the same scenario to evaluate the performance of our NNCG solver, in comparison to a default Jacobi solver as used by Gissler et al. [GPB*19]. For that, we simulate 1 s of the blocks dropping and resting on each other. In both cases, we set t = 0.001s, υ = 0.5 with particle distance h = 0.02 m and the same desired residual. As we can see in the top diagram in Figure 7, the default Jacobi solver on average requires approximately five times more iterations compared to the NNCG solver. This speedup by far overcompensates for the slightly higher computational cost per NNCG iteration, resulting in a significant Figure 5: Degenerated test cases proposed by Erleben [Erl18] which can be challenging for mesh-based contact handling methods. For each of the nine constellation, we show the initial setting and right next to it a snapshot of the simulation. The bodies are given some initial velocity as described in Erleben [Erl18]. Visually, our particle-based contact handling shows no problems whatsoever when simulating the shown constellations. Figure 6: A heavy rigid block (yellow) is dropped on a much lighter block (red) where it then comes to rest. The mass ratio between the blocks is 1000 : 1. Our contact handling method is able to robustly handle this challenging constellation.
performance advantage for the NNCG solver in this scenario. However, we note that this scenario represents a special case of a very challenging body constellation, which is why we additionally evaluate the solver performance at a more common and better behaved body arrangement. We simulate a stack of eight rigid boxes with equal masses using the same parameter setting as before, and again compare the solver iterations of the NNCG solver against the iterations required by the Jacobi solver. The result is visualized in the centre of Figure 7. As we can see, the difference in the number of solver iterations is a little less severe compared to the previous scenario and in some simulation steps, the NNCG solver takes almost as many iterations as the Jacobi solver. However, conforming to the previous scenario, on average, the Jacobi solver still requires more than three times as many iterations as the NNCG solver when simulating the stacked boxes. To evaluate the solver performance in a coupled fluid-rigid simulation, we simulate a fluid pillar with a rigid duck swimming on the surface. State-of-the-art particle fluid solvers such as PCISPH [SP09], PBF [MM13], IISPH [ICS*14] and DF-SPH [BK15] use a Jacobi-style solving method comparable to the In the example with the high body mass ratio, at the initial contact, the Jacobi solver requires approximately 10 times more iterations compared to the NNCG solver to achieve the same residual error. After the bodies came to rest, the Jacobi solver on average still requires 5 times as many iterations. To simulate the stacked boxes, the Jacobi solver requires a little less than 4 times as many solver iterations, and in the fluid pillar example, the Jacobi solver even requires over 10 times more iterations. In all three test scenarios, our NNCG solver has a significant performance advantage over a standard Jacobi method.
one presented by Gissler et al. [GPB*19]. We evaluate our implementation in relation to a standard Jacobi solver similar to Gissler et al. [GPB*19] in order to classify the performance of the NNCG solver in the context of fluid and coupled fluid-rigid simulations. The required solver iterations of the Jacobi and NNCG solver are displayed at the bottom of Figure 7. The difference in number of solver iterations in this scenario is especially significant as the Jacobi solver requires more than 10 times the iterations the NNCG solver needs to reach an equal residual.
Motivated by the presented performance advantages, as a last test scenario, we evaluate how the Jacobi and NNCG solver scale for increasing scenario complexity. A commonly used factor to determine scenario complexity is depth of a fluid body under gravity [KBST19]. We simulate a fluid pillar that grows in height during the simulation and measure the required solver iterations over time. The results are displayed in Figure 8. As we can see, the Jacobi iterations increase approximately linearly with increasing fluid depth, while the number of NNCG iterations show sub-linear growth. Thus, in this specific scenario, for rising fluid depth, the significance of the performance advantage of the NNCG solver increases. In summary, in all tested scenarios, the NNCG solver maintains a considerable performance advantage over the Jacobi solver.

Versatility and scalability
We demonstrate the capabilities of our simulation method by showcasing a number of simulation scenarios we deemed to be challenging to the friction and pressure solver due to many interacting rigid bodies, complex interactions between fluids and rigid geometry as well as scenarios requiring precisely computed friction forces. For each scenario, average timesteps, required solver iterations and computation times are listed in Table 1.

Lewis lifting mechanism: Inspired by Ferguson et al. [FLS*21]
, we simulate a Lewis lift which is a mechanism that can lift weights relying on frictional forces. As Figure 9 shows, our friction solver is able to handle static friction in a correct and precise manner without revealing any visible artifacts.
Modified card house: Simulating houses of cards is a popular way to demonstrate correctness and precision of frictional forces [KSJP08, BET14, LFS*20, FLS*21]. To show that our friction handling matches the ability of state-of-the-art friction computation methods, we simulate a seven-story card house. To further increase difficulty, we added a hole in the middle of the card house making the structure even more fragile. Figure 10 shows the card house and its destruction.
Breaking dam: A dam consisting of stacked rigid blocks is holding back a body of water whose depth is continuously rising. Here, the strong coupling between fluid pressure forces and dry friction is especially important to ensure a stable and robust simulation at the It is easy to see that our proposed strong coupling between fluid pressure and rigid friction forces helps to guarantee a stable simulation of the water-dam interface.
interface between fluid and dam. As soon as pressure forces caused by the fluid are growing larger than the static friction forces between the blocks can compensate for, the dam starts to break. As illustrated in Figure 11, the released fluid causes a flood and parts of the dam knock over a bridge standing downstream. During the simulation, the fluid is sampled by a maximum of 2.5 million particles, the dam and bridge are build from 128 rigid blocks sampled by 600 thousand particles and the kinematic boundary is sampled with 1.2 million particles.
Tower: Simulating high stacks of rigid bodies is inherently difficult as contact information needs to be propagated through the stack to ensure incompressibility of all rigid bodies and correct friction handling [Erl07]. Our solver is able to stably simulate an over 40story high tower standing upright. By shooting two rigid bodies into the tower, we let it collapse into a basin filled with fluid, showcasing the robustness of the monolithic simulation method. Figure 12 gives an overview of the simulation using 3.6 million fluid particles, Figure 12: Over 600 rigid blocks are used to build a tower. Frictional forces ensure that the tower is able to stand upright, even with fluid pushing against the base of the tower. Using a waterslide, two rigid ducks crash into the tower causing it to collapse.
1.4 million rigid particles sampling 658 rigid bodies and 4.9 million particles belonging to kinematic boundaries.

Bulk simulation:
We showcase a bulk simulation to hint at the wide range of possible applications of our simulation method. As shown in Figure 13, the bulk consists of thousands of rigid bodies in the shape of armadillos. Each armadillo is sampled with 622 rigid particles, giving a good approximation of the complex underlying geometry. To demonstrate the robustness of our solver when confronted with complex interacting geometry under difficult conditions, the armadillos are caught in a net consisting of actively simulated rigid chain links. Over time, a pile of armadillos is formed which is held in place by static friction forces. In total, there are 8700 armadillos and 670 chain links sampled with 5.7 million rigid particles. The kinematic geometry is sampled using 5.0 million particles.
Chain: Simulating long chains of rigid bodies is inherently difficult due to mutual dependencies of individual chain links over long distances. Contacts need to be resolved very precisely to prevent elongations and elastic behaviour of the chain. Illustrated in Figure 14, we demonstrate that our solver is able to simulate a 100 rigid links long chain, with no need for an extra long-range constraint stabilization technique such as proposed in Müller et al. [MCMJ17]. Additionally, to further increase difficulty, some chain links are shaped as hollow letters containing fluid and small rigid ducks. The chain is rolled up on an axle causing challenging simulation circumstances with high pressure and friction forces between contacting chain links, diverse sizes and geometries of links and rapid velocity changes in the chain including the contained fluid. In total, the chain consists of over 100 rigid bodies, sampled by 600 thousand particles. The fluids are represented using 270 thousand particles.

Limitations
Since the number of particles used to sample rigid body geometry mostly depends on body and particle size instead of complexity of the geometry, simple shapes with large flat surfaces still might require thousands of sample particles. In this situation, a pure rigid body simulation ignoring rigid-fluid interaction is likely to be more efficient when using a mesh-based rigid body representation. Our simulation method uses the same size for rigid body particles, fluid particles and boundary particles. While this allows a straightforward interaction between particles, the resolution of the contact handling between rigid bodies is directly coupled to the fluid resolutions. In the future, approaches that implement different particle sizes such as described in Refs. [WHK17, WK21] could be investigated to enable a more flexible simulation resolution. Due to the velocity-based friction formulation drift might occur between contacting rigid bodies even if static friction forces are applied. This drift can be reduced at the cost of more solver iterations, but it cannot be eliminated entirely. In the future, it might be interesting to additionally implement a stabilization scheme that guarantees static contacts such as described by Xu et al. [XZB14] or Ding and Schroeder [DS20]. In our simulation, we employed a simple Coulomb friction model. However, we believe that in the future, the described solver could be easily extended to handle more sophisticated friction models such as anisotropic friction [Erl17,EMAK19] or torque due to dry friction resisting spinning [BNT*15].

Conclusion
We have introduced a monolithic SPH solver for particle-based fluids and rigid bodies including dry frictional forces. Our simulation method calculates predicted volume errors as a simple and robust way to detect collisions between rigid bodies as well as compressed fluid particles. Volume errors are resolved by computing pressure and pressure forces, a concept well known from existing particlebased pressure solvers [SP09, ICS*14, BK15, GPB*19]. Strong coupling between pressure and friction forces is achieved by directly incorporating the effects of pressure and fictional forces into the computation of predicted velocities, which in turn influence the current pressure and friction computation. A friction mirroring step is build into the pressure and friction system to guarantee conservation of momentum. Building on top of the relaxed-Jacobi solver, a non-linear non-smooth conjugate gradient method is employed to accelerate the solving procedure. Depending on the scenario, this can result in a significant reduction of solver iterations. Further, we have shown that our implicit friction formulation can handle static friction, produces correct behaviour at stick-slip transitions and scales to large numbers of simultaneous rigid body contacts. The strong coupling between pressure and friction allows us to stably simulate complex structures of rigid bodies interacting with fluids, opening the door for a whole range of new interesting simulation scenarios.
Algorithm 6. The procedure to efficiently compute the relation between pressure p r and volume error V err r (t + t) for a rigid particle r. All quantities that are explicitly computed in the algorithm are pre-computed and used again later while quantities that are not mentioned are computed on-the-fly.
The term ∂/∂ λ r v rel r (t + t ) in Equation (A.10) is also constant during one simulation step. Using Equation (21a), we can write where velocities can be written as described in Equation (11): Note that we switched the order of the cross products to express them as matrix multiplications. Frictional forces F F and torque τ F are computed in the friction mirroring step described in Equation (22). This gives us translational velocity v R and angular velocity ω R using We use quaternions to represent orientations and employed a time integration scheme proposed by Zhao and Wachem [ZW13] that is able to integrate quaternions with no need for renormalization. Afterwards, we compute the new particle positions and velocities directly based on x R (t + t ), q R (t + t ), v R (t + t ) and ω R (t + t ): x r (t + t ) = x R (t + t ) Algorithm 9. The implementation of the Euler-Cromer time integration scheme used in our simulation.
v r (t + t ) = v R (t + t ) Algorithm 9 summarizes the time integration step.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.