Strongly coupled simulation of incompressible fluid and rigid bodies with velocity-based constraints using particle method

TL;DR

This paper introduces a particle-based, velocity-formulated method for strongly coupling incompressible fluid and rigid bodies, ensuring compatibility of constraints and enabling complex simulations.

Contribution

It presents a novel velocity-based formulation and a unified approach to simulate incompressible fluids and rigid bodies with strong coupling using particle methods.

Findings

Successfully simulates dam-break and hydrostatic pressure scenarios.

Demonstrates stable interaction between fluid and rigid bodies.

Validates the method with complex scene simulations.

Abstract

This paper presents a novel particle method to compute strongly coupled incompressible fluid and rigid bodies. The method adopts a velocity-based formulation and utilizes the linear complementarity problem for the incompressibility constraint. Since all the constraints for incompressibility, inter-rigid-body contacts, and interaction between incompressible fluid and rigid bodies are mathematically compatible, strongly coupled simulation is achieved using the method, where the shapes of the rigid bodies are represented by particles as well. The abstract concept of velocity-based constraints is presented, which generalizes the formulations of the incompressibility constraint and inter-rigid-body contacts and provides a generic way to achieve strongly coupled simulation. Several numerical examples are presented to verify the method, which includes rigid-body computation, hydrostatic pressure, dam-break computation, and circular parch computation for incompressible fluid, buoyancy and seesaw computation for interaction of incompressible fluid and rigid bodies, and complex-scene computation for overall behavior and stability.

Equations

$$\frac{D\mathbf{u}}{Dt}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}+\mathbf{f}$$

$$\nabla \cdot \mathbf{u}=0$$

$$w(r)=\{\begin{array}{ll}(1-\frac{r}{r_e}{)}^2& (r<r_e)\\ 0& (\text{otherwise})\end{array}$$

$$n_i=\sum _{j\ne i}w(r_{ij}),$$

$${⟨\nabla \varphi ⟩}_i=C\sum _{j\ne i}({\varphi }_i+{\varphi }_j)w^{\prime }(r_{ij})\frac{{\mathbf{r}}_{ij}}{r_{ij}},$$

$$C={\left(\sum _{j\ne i}({\mathbf{r}}_i+{\mathbf{r}}_j{)}_x{w}^{\prime }(r_{ij})\frac{({\mathbf{r}}_{ij}{)}_x}{r_{ij}}\right)}^{-1}$$

$$-\frac{1}{\rho }⟨\nabla p{⟩}_i=-\frac{C}{\rho }\sum _{j\ne i}(p_i+p_j)w^{\prime }(r_{ij})\frac{{\mathbf{r}}_{ij}}{r_{ij}}.$$

$$n_i=n^0.$$

$$\frac{d{n}_i}{dt}=0.$$

$$\frac{d{n}_i}{dt}$$

$$=\sum _{j\ne i}{w}^{\prime }(r_{ij})\frac{d}{dt}\Vert {\mathbf{r}}_j-{\mathbf{r}}_i\Vert$$

$$=\sum _{j\ne i}{w}^{\prime }(r_{ij})\frac{({\mathbf{r}}_j-{\mathbf{r}}_i)\cdot ({\mathbf{u}}_j-{\mathbf{u}}_i)}{\Vert {\mathbf{r}}_j-{\mathbf{r}}_i\Vert }$$

$$=\sum _{j\ne i}{w}^{\prime }(r_{ij})({\mathbf{u}}_j-{\mathbf{u}}_i)\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}},$$

$$\sum _{j\ne i}{w}^{\prime }(r_{ij})({\mathbf{u}}_j-{\mathbf{u}}_i)\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}}=0.$$

$$\sum _{j\ne i}{w}^{\prime }(r_{ij})({\mathbf{u}}_j-{\mathbf{u}}_i)\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}}$$

$$\frac{\alpha }{h}(n^0-n_i)$$

$$=\sum _{j\ne i}{w}^{\prime }(r_{ij})\left({\mathbf{u}}_j^{\ast }-{\mathbf{u}}_i^{\ast }-\frac{h}{\rho }⟨p_j^{t+1}⟩+\frac{h}{\rho }⟨p_i^{t+1}⟩\right)\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}},$$

$$\frac{\alpha }{h}(n^0-n_i)-\sum _{j\ne i}{w}^{\prime }(r_{ij})({\mathbf{u}}_j^{\ast }-{\mathbf{u}}_i^{\ast })\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}}$$

$$=h\frac{C}{\rho }\sum _{j\ne i}{w}^{\prime }(r_{ij})\left(\sum _{k\ne i}(p_i+p_k)w^{\prime }(r_{ik})\frac{{\mathbf{r}}_{ik}}{r_{ik}}-\sum _{k\ne j}(p_j+p_k)w^{\prime }(r_{jk})\frac{{\mathbf{r}}_{jk}}{r_{jk}}\right)\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}}.$$

$${\mathbf{v}}_{ij}=\{\begin{array}{ll}0& (i=j)\\ w^{\prime }(r_{ij})\frac{{\mathbf{r}}_{ij}}{r_{ij}}& (i\ne j).\end{array}$$

$$h\frac{C}{\rho }\sum _{j\ne i}\sum _{k\ne i}(p_i+p_k)w^{\prime }(r_{ij})\frac{{\mathbf{r}}_{ij}}{r_{ij}}\cdot w^{\prime }(r_{ik})\frac{{\mathbf{r}}_{ik}}{r_{ik}}$$

$$-h\frac{C}{\rho }\sum _{j\ne i}\sum _{k\ne j}(p_j+p_k)w^{\prime }(r_{ij})\frac{{\mathbf{r}}_{ij}}{r_{ij}}\cdot w^{\prime }(r_{jk})\frac{{\mathbf{r}}_{jk}}{r_{jk}}$$

$$=h\frac{C}{\rho }\sum _{j,k}(p_i+p_k){\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{ik}-h\frac{C}{\rho }\sum _{j,k}(p_j+p_k){\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{jk}$$

$$\displaystyle=h\frac{C}{\rho}\biggl{(}p_{i}\sum_{j,k}{\boldsymbol{v}}_{ij}\cdot{\boldsymbol{v}}_{ik}+\sum_{j}p_{j}\sum_{k}{\boldsymbol{v}}_{ij}\cdot{\boldsymbol{v}}_{ik}$$

$$\displaystyle\phantom{=}-\sum_{j}p_{j}\sum_{k}{\boldsymbol{v}}_{ij}\cdot{\boldsymbol{v}}_{jk}+\sum_{j}p_{j}\sum_{k}{\boldsymbol{v}}_{ik}\cdot{\boldsymbol{v}}_{jk}\biggr{)}$$

$$=p_{i}h\frac{C}{\rho }\left({\Vert \sum _k{\mathbf{v}}_{ik}\Vert }^2+\sum _k{\Vert {\mathbf{v}}_{ik}\Vert }^2\right)$$

$$+\sum _{j\ne i}{p}_{j}h\frac{C}{\rho }\sum _k\left({\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{ik}-{\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{jk}+{\mathbf{v}}_{ik}\cdot {\mathbf{v}}_{jk}\right).$$

$$A_{ij}=h\frac{C}{\rho }\{\begin{array}{ll}{\Vert \sum _k{\mathbf{v}}_{ik}\Vert }^2+\sum _k{\Vert {\mathbf{v}}_{ik}\Vert }^2& (i=j)\\ \sum _k\left({\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{ik}-{\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{jk}+{\mathbf{v}}_{ik}\cdot {\mathbf{v}}_{jk}\right)& (i\ne j).\end{array}$$

$${\mathbf{J}}_i=\left(\begin{array}{c}{\mathbf{v}}_{i,1}\\ ⋮\\ {\mathbf{v}}_{i,i-1}\\ -{\sum }_j{\mathbf{v}}_{i,j}\\ {\mathbf{v}}_{i,i+1}\\ ⋮\\ {\mathbf{v}}_{i,N}\end{array}\right),$$

$$h\frac{C}{\rho }{\mathbf{J}}_i\cdot {\mathbf{J}}_j$$

$$=h\frac{C}{\rho }\{\begin{array}{ll}\sum _k{\Vert {\mathbf{v}}_{ik}\Vert }^2+{\Vert \sum _k{\mathbf{v}}_{ik}\Vert }^2& (i=j)\\ \sum _k({\mathbf{v}}_{ik}\cdot {\mathbf{v}}_{jk}-{\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{jk}+{\mathbf{v}}_{ij}\cdot {\mathbf{v}}_{ik})& (i\ne j)\end{array}$$

$$=A_{ij},$$

$$(\mathbf{p}{)}_i$$

$$(\mathbf{A}\mathbf{p}-\mathbf{b}{)}_i$$

$$(\mathbf{p}{)}_i(\mathbf{A}\mathbf{p}-\mathbf{b}{)}_i$$

$${\delta }_i=\frac{\alpha }{h}(n^0-n_i)-\sum _{j\ne i}{w}^{\prime }(r_{ij})({\mathbf{u}}_j^{\ast }-{\mathbf{u}}_i^{\ast })\cdot \frac{{\mathbf{r}}_{ij}}{r_{ij}}.$$

$$p_i^{\prime }=\max \left\{0,p_i^{t+1}+\frac{{\delta }_i}{A_{ii}}\right\}.$$

$$\Delta p_i$$

Full text

1 Introduction

Fluid simulation using computers has widely been used for various purposes, such as analysis of industrial products, generating computer graphics, and usage in computer games. Although the dynamics of fluid can be described by the Navier-Stokes equations, solving these equations analytically is unrealistic in many cases and, therefore, solving these equations numerically using computers is important. In reality, fluid interacts with various solids. As such, enabling fluid-solid interaction makes simulation more practical. In addition, within the range of daily handling, many types of fluids and solids can be approximated as incompressible fluids or rigid bodies, respectively. Therefore, coupling of incompressible fluid and rigid body enables more versatile and complex simulation. Coupled simulation of multiple substances can be widely divided into two types of simulation, namely, weakly coupled simulation and strongly coupled simulation. In weakly coupled simulation, interaction between different substances is computed explicitly, whereas, in strongly coupled simulation, this interaction is implicitly computed. Therefore, strongly coupled simulation is more stable and accurate than weakly coupled simulation when the same time step size is adopted.

Rigid-body simulation has long been researched, and impulse-based methods [1] are among the main computation methods. The most important part of rigid-body computation is contact computation, and impulse-based methods deal with contacts in a unified manner by the impulse that integrates force over a small time step. Since simulation requires the integration of time, impulse is easier to handle than force, and using impulse keeps simulation from mathematical and numerical problems [2, 3]. If we only consider normal force, we can formulate contact constraints using linear complementarity problems (LCPs) [4]. When we take friction into consideration, the problem is no longer linear and we have to solve nonlinear complementarity problems (NCPs). However, we can approximate friction to make things less complicated. For example, Tonge et al. [5] formulated frictional constraints using boxed LCPs by applying pyramid approximation to the Coulomb friction cone. Moreover, Gholami et al. [6] adopted a continuity approximation to friction to formulate the entire contact problem using LCPs. On the other hand, we can also compute the accurate Coulomb friction cone without approximation by applying an extension to the projected Gauss-Seidel (PGS) method, which is generally used to solve LCPs [7]. We also use the PGS method and its extension for solving all kinds of constraints that appear in simulation.

When we formulate rigid-body contacts using LCPs, we have to predict the relative velocity after collisions prior to contact computation. This prediction is usually performed using the relationship between the relative velocity and the coefficient of restitution. However, when multiple collisions occur simultaneously, it is difficult to explicitly compute the accurate velocity after these collisions, and LCPs may give physically incorrect results. In addition, depending on the positions and number of contact points, the coefficient matrix of an LCP may not become a symmetric positive-definite matrix, and therefore the solution of an iterative method may not converge to a unique solution [8]. Although we do not adopt them in this study, some approaches are presented to solve these problems. For example, Tang et al. [9] developed the energy tracking impulse (ETI) method, which does not require solving LCPs to compute the impulse. In the ETI method, by tracking energy during collisions, the physically accurate impulse can be calculated without explicitly giving the relative velocity after collision, and its improvement was proposed by Li et al. [10] for particle-based rigid-body simulation.

Linear complementarity problems have long been used for rigid-body computation, but there are a few examples of the use of LCPs for fluid computation. Batty et al. [11] used an LCP to formulate wall boundary conditions, and Bodin et al. [12] proposed a method by which to simultaneously formulate and solve constant-density constraints and constraints for wall boundary conditions as a mixed LCP. In addition, Gerszewski and Bargteil [13] developed a method by describing the incompressibility of a fluid using inequalities and formulating the resulting pressure equation with an LCP.

Coupled simulation of fluid using a particle method and rigid bodies has long been researched. In the context of the smoothed particle hydrodynamics (SPH) method, Monaghan et al. [14] developed a method that enables coupled simulation of fluid and rigid bodies by representing rigid bodies with particles and computing the interaction force between nearing fluid particles and rigid-body particles. Later, a more stable method was developed by Oh et al. [15] to compute the interaction based on impulse rather than force, including inter-rigid-body collisions. Akinci et al. [16] improved computation at the boundaries of a particle-based rigid-body and fluid and generalized the explicit computation of the interaction so that they can handle the coupled simulation of rigid bodies with fluid computed by both the predictive-corrective SPH method [17] and the slightly compressible SPH method. Macklin et al. [18] proposed a method that enables coupled simulation of liquid, gas, solid, and cloth on GPUs by approximating rigid bodies using shape matching techniques. Coupling of rigid bodies and an incompressible fluid has also been researched. Li and Asai [19] proposed a method to couple an incompressible fluid using the incompressible SPH method [20] and rigid-bodies that are computed using an impulse-based method. There is also a method, which was developed by Klinger et al. [21], to strongly couple an incompressible fluid and rigid bodies using meshes that are dynamically updated during simulation. In the field of computer graphics, incompressible fluid simulation with position-based formulation has been proposed [22], which can easily be coupled with other position-based methods. However, position-based methods cannot capture some types of physical phenomena, such as variation of coefficient of restitution in rigid-body simulation.

In addition, the moving particle semi-implicit (MPS) method [23] was used to research the coupling of fluids and rigid bodies. Shibata et al. [24] simulated shipping water on a ship to analyze its effect using weak coupling of fluid and a rigid body. In addition, Shibata et al. [25] simulated a lifeboat falling into water also using weak coupling of fluid and a rigid body by transferring fluid pressure to the rigid body. Koshizuka et al. [26] developed the passively moving solid (PMS) model, which computes rigid bodies with particles by once treating rigid-body particles as fluid particles during pressure and advection computation and then cancelling out the deformation of rigid bodies to restore their shape. The PMS model can handle not only the interaction between rigid bodies and fluid but also between rigid bodies and other rigid bodies. Tanaka et al. [27] computed coupled simulation of rigid bodies and an incompressible fluid using a penalty method for inter-rigid-body contact calculation and the PMS model for the interaction of rigid bodies and fluid. Gotoh et al. [28] computed a flood flow with floating objects using the MPS method and the PMS model for fluid-and-solid coupled simulation. The PMS model itself is not necessarily combined with the MPS method, and the same approach has been used with the SPH method [29] and the weakly compressed SPH method [30] to compute floating rigid bodies. This approach has also been used in the field of computer graphics [31]. Large-scale computation of fluid and rigid bodies has been performed by Murotani et al. [32] using the explicit MPS method for fluid and the PMS model for the interaction of fluid and rigid bodies. They divided the computation space into multiple small spaces in order to make the simulation run in parallel [32]. Wang et al. [33] simulated floating bodies transported by fluid using the MPS method and the PMS model for an incompressible fluid and the interaction of rigid bodies and fluid, as well as the discrete element method for explicitly calculated inter-rigid-body contacts. In the MPS method, the pressure of a fluid is implicitly computed by solving the pressure Poisson equation and thus the pressure field of rigid bodies in the PMS method is implicitly calculated along with that of the fluid when this method is used in combination with the MPS method. As far as we know, there is no other method than ours that achieves fully strongly coupled simulation of incompressible fluid and rigid bodies using a particle method and velocity-based formulations.

Rigid-body computation in impulse-based methods is performed in a velocity-based manner. Therefore, by formulating the incompressibility constraint of fluid in a velocity-based manner as well, we can expect to simultaneously compute both the impulse due to the collision of rigid bodies and the pressure of an incompressible fluid implicitly [34]. In this study, we modify and generalize our previous work [34]. Although the strongly coupled simulation itself has already been achieved in [34], the core technique of the method highly depends on LCP formulations. In the present paper, we extend the core technique that enables the strongly coupled simulation by introducing the abstract concept of “velocity-based constraints” so that it can be applied for wider range. The notion of velocity-based constraints provides a very flexible framework of strongly coupled simulation of various substances that are not limited to fluids and rigid bodies.

The remainder of the present paper is organized as follows. In Section 2, we introduce the proposed particle method to compute incompressible fluid. We also define the velocity-based constraints in this section. In Section 3, we describe the impulse-based rigid-body computation and the way to build velocity-based constraints to solve inter-rigid-body collisions. Then, in Section 4, the strong coupling of rigid bodies and an incompressible fluid is proposed. We present some numerical examples in Section 5 in order to confirm the behavior and accuracy of the proposed method. Finally, we present the conclusion in Section 6.

2 Fluid Simulation

In this section, we introduce a velocity-based method to simulate incompressible fluids. We use ${\boldsymbol{r}}$, ${\boldsymbol{u}}$, $m$, $\rho$, $p$, $\nu$, and ${\boldsymbol{f}}$ to represent position, velocity, mass, density, pressure, kinematic viscosity, and acceleration due to external forces, respectively. We use subscript $\phi_{i}$ to represent arbitrary physical quantity $\phi$ of particle $i$.

2.1 Governing Equation

We use the Navier-Stokes equation

[TABLE]

and the continuity equation

[TABLE]

as the governing equations of incompressible fluids. Here, we assumed that the density and the viscosity are constant. Note that $D{\boldsymbol{u}}/Dt$ in the left-hand side of (1) represents the Lagrange derivative of ${\boldsymbol{u}}$.

2.2 Weighting Function

We denote the weighting function as $w$ and the effective radius as $r_{e}$. In the proposed method, the derivative of weighting function $w^{\prime}$ is used during computation, so it is preferable that $w^{\prime}$ be continuous around effective radius $r_{e}$. Thus, we use the following function as the weighting function:

[TABLE]

Using the weighting function, we define the particle number density of particle $i$ as

[TABLE]

where $r_{ij}={\left\lVert{\boldsymbol{r}}_{ij}\right\rVert}$ and ${\boldsymbol{r}}_{ij}={\boldsymbol{r}}_{j}-{\boldsymbol{r}}_{i}$. We also define the standard particle number density $n^{0}$ as the particle number density in the initial particle arrangement of a particle that is placed inside a sufficient volume of the fluid. We use $l$ to denote the interval between particles in the initial particle arrangement.

2.3 Spatial Discretization

We define the discretized gradient model as

[TABLE]

which is widely used in various SPH-based methods, except for the scalar constant $C$. We introduce the constant because the most standard gradient model used in the SPH method tends to give an incorrect value, especially when the effective radius is not large enough. The value of $C$ is chosen so that the gradient model shows the correct value in the initial particle arrangement. More precisely, $C$ is computed as

[TABLE]

in the initial particle arrangement, where $(\cdot)_{x}\colon{\mathbb{R}^{3}}\to\mathbb{R}$ is the first component of a vector and particle $i$ is inside a sufficient volume of the fluid.

Using the gradient model of (5), the pressure term for particle $i$ is discretized as follows:

[TABLE]

2.4 Incompressibility Constraint

Since the particle number density is proportional to the fluid density, we can adopt the following equation as the incompressibility condition:

[TABLE]

From (8), we obtain

[TABLE]

By the definition of the particle number density, we can directly calculate the left-hand side of (9) as:

[TABLE]

and thus the discretized version of the zero-divergence condition of (2) can be written as:

[TABLE]

However, solving only (14) as the constraint in actual computation makes the numerical error accumulate as the computation progresses, which allows the fluid to be gradually compressed, and eventually the computation collapses. There are several approaches to avoid this problem. One approach is to utilize the particle shifting techniques [35, 36] to maintain the volume of the fluid constant. Another approach is to mix the positional constraint violation term to the velocity constraint. The latter is easier to handle because this approach requires few changes to the entire procedure and keeps the equations to be solved linear. We use this approach to solve the problem. Let $\alpha\in[0,1]$ be a constant to control the amount of mixing. We use the following equation to obtain the constraint:

[TABLE]

where $h$ is the time step size. Setting $\alpha$ to zero means that there is no positional term and (15) is equal to (14), whereas setting $\alpha$ to one means that an attempt is made to cancel out all positional constraint violations at the next time step. However, using a large value of $\alpha$ does not work well due to the linear approximation and makes the simulation unstable, so we use 0.05 for the value of $\alpha$ throughout the present paper.

2.5 Time Integration

In this method, we use a semi-implicit scheme for the time integration, which is similar to the MPS method. The computation in each time step is performed roughly as follows. We first calculate the viscosity and the external force terms explicitly to obtain temporal velocity ${\boldsymbol{u}}^{*}$, and then we compute the pressure to correct the temporal velocity so that ${\boldsymbol{u}}^{*}$ satisfies constraint  eq:fluid:mixed. Finally, we adopt the corrected temporal velocity as the velocity in the next time step, and update the position using the velocity.

We use integer $t$ to denote a time step, and use superscript $\phi^{t}$ to denote arbitrary physical quantity $\phi$ at time $t$. We do not define the detailed algorithm for computing explicit terms, as there are many acceptable ways to compute these terms and the subsequent procedure is not affected as long as we can obtain temporal velocity ${\boldsymbol{u}}^{*}$. In the following discussion, we focus on the computation of ${\boldsymbol{u}}^{t+1}$ and $p^{t+1}$.

Let ${\boldsymbol{x}}^{t}$ and ${\boldsymbol{u}}^{t}$ be the position and the velocity of particles, respectively, at time $t$. After obtaining temporal velocity ${\boldsymbol{u}}^{*}$ by explicit calculation, we first initialize pressure $p^{t+1}$ with zero. Pressure $p^{t+1}$ plays a role in correcting temporal velocity ${\boldsymbol{u}}^{*}$ through the pressure term $-\frac{1}{\rho}\left\langle p^{t+1}\right\rangle$. In order to derive an equation to obtain $p^{t+1}$, assume that ${\boldsymbol{u}}^{t+1}={\boldsymbol{u}}^{*}-\frac{h}{\rho}\left\langle p^{t+1}\right\rangle$ holds and assign ${\boldsymbol{u}}^{t+1}$ to constraint (15) to obtain

[TABLE]

or, equivalently,

[TABLE]

Let $N$ be the number of particles; ${\boldsymbol{b}}\in{\mathbb{R}^{N}}$ be a vector, the $i$th component of which is $\frac{\alpha}{h}(n^{0}-n_{i})-\sum_{j\neq i}w^{\prime}(r_{ij})({\boldsymbol{u}}_{j}^{*}-{\boldsymbol{u}}_{i}^{*})\cdot\frac{{\boldsymbol{r}}_{ij}}{r_{ij}}$; and ${\boldsymbol{p}}\in{\mathbb{R}^{N}}$ be a vector, the $i$th component of which is $p_{i}^{n+1}$. Then, there exists a matrix ${\boldsymbol{A}}\in{\mathbb{R}^{N\times N}}$ such that ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ is equivalent to equation (19). Define ${\boldsymbol{v}}_{ij}\in{\mathbb{R}^{3}}$ as

[TABLE]

Then, the right-hand side of (19) can be written as

[TABLE]

Therefore, the components of matrix ${\boldsymbol{A}}$ are

[TABLE]

Let ${\boldsymbol{J}}_{i}\in{\mathbb{R}^{3N}}$ be a vector

[TABLE]

and ${\boldsymbol{J}}\in{\mathbb{R}^{3N\times N}}$ be a matrix ${\boldsymbol{J}}=({\boldsymbol{J}}_{1}\cdots{\boldsymbol{J}}_{N})$. Then, matrix ${\boldsymbol{A}}$ can be decomposed into ${\boldsymbol{A}}=h\frac{C}{\rho}{\boldsymbol{J}}^{T}{\boldsymbol{J}}$. Indeed, the inner product of ${\boldsymbol{J}}_{i}$ and ${\boldsymbol{J}}_{j}$ times $h\frac{C}{\rho}$ is

[TABLE]

and, therefore, if the particles are well arranged, in other words, if matrix ${\boldsymbol{J}}$ is not degenerated, then matrix ${\boldsymbol{A}}$ is symmetric positive definite.

As discussed above, we need to solve the equation ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ with respect to ${\boldsymbol{p}}$ in order to correct the temporal velocity so that constraint (15) is satisfied. However, directly using the pressure obtained by solving the equation may make the simulation unstable. Particle methods sometimes behave poorly when there is an attractive force between particles. Considering the gradient model of (5), we see that the attractive force due to the pressure works between particles $i$ and $j$ if and only if $p_{i}+p_{j}<0$ holds for some $i$ and $j$. This fact leads us to an easy sufficient condition: there is no attractive force if $p_{i}\geq 0$ holds for every particle $i$. The simplest way to make the pressure nonnegative is just to project the solution of ${\boldsymbol{A}}{\boldsymbol{x}}={\boldsymbol{b}}$ onto $\mathbb{R}_{+}^{N}$, where $\mathbb{R}_{+}=[0,\infty)$. However, this modification of the solution breaks the equilibrium between the pressure of each particle. If the pressure of particle $i$ is changed from negative to zero, then the $j$th row of equation ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ will no longer hold for each particle $j$ that is a neighboring particle of particle $i$. We need ${\boldsymbol{p}}$ to satisfy ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ as much as possible, whereas ${\boldsymbol{p}}$ is forced to be nonnegative. These requirements can be satisfied by finding ${\boldsymbol{p}}$ that satisfies

[TABLE]

for each $i$, where $(\cdot)_{i}$ is the $i$th component of a vector. Inequality (33) is for the nonnegativity of ${\boldsymbol{p}}$. Inequality (34) is for the requirement that ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ should hold as much as possible. If $({\boldsymbol{A}}{\boldsymbol{p}}-{\boldsymbol{b}})_{i}<0$ holds, then we can make $({\boldsymbol{A}}{\boldsymbol{p}}-{\boldsymbol{b}})_{i}=0$ hold by increasing $({\boldsymbol{p}})_{i}$, because ${\boldsymbol{A}}$ is positive definite and, therefore, $A_{ii}>0$. Equation (35) has the same requirement that ${\boldsymbol{A}}{\boldsymbol{p}}={\boldsymbol{b}}$ should hold as much as possible. If $({\boldsymbol{p}})_{i}({\boldsymbol{A}}{\boldsymbol{p}}-{\boldsymbol{b}})_{i}>0$ holds, then $({\boldsymbol{p}})_{i}$ is unnecessarily large, and either $({\boldsymbol{p}})_{i}=0$ or $({\boldsymbol{A}}{\boldsymbol{p}}-{\boldsymbol{b}})_{i}=0$ can be achieved by decreasing $({\boldsymbol{p}})_{i}$ because $A_{ii}>0$ holds. The problem of finding ${\boldsymbol{p}}$ that satisfies (33), (34), and (35) is called a linear complementarity problem (LCP). An LCP has a unique solution if the coefficient matrix of the problem is symmetric positive definite. Since ${\boldsymbol{A}}$ is proven to be symmetric positive definite, there exists a unique solution of the problem.

We configured matrix ${\boldsymbol{A}}$ in order to define the problem. However, in actual computation, we do not compute the matrix explicitly. Instead, we use temporal velocity ${\boldsymbol{u}}^{*}$ to iteratively solve this model. Since the constraint is originally imposed on the velocity, we repeatedly apply the change in pressure to the temporal velocity using the gradient model of (7), while determining the change in pressure by computing the current deviation of ${\boldsymbol{u}}^{*}$ from constraint (15). More precisely, we iterate the following steps until some stopping criteria are satisfied. For each particle $i$, we first compute the current deviation of constraint $\delta_{i}$ as

[TABLE]

Then, we compute the next value of pressure $p^{\prime}_{i}$ as

[TABLE]

Finally, we update pressure $p_{i}^{t+1}$ and apply the difference of it to temporal velocity ${\boldsymbol{u}}^{*}$:

[TABLE]

This iterative approach actually corresponds to the projected Gauss-Seidel (PGS) method, which is an iterative method to solve LCPs. The tentative solution is updated repeatedly in a Gauss-Seidel manner and is projected onto the nonnegative space.

Remarkably, we can compute pressure $p^{t+1}$ based only on temporal velocity ${\boldsymbol{u}}^{*}$. When we update the pressure of a particle, we no longer need to refer the pressure of the neighboring particles because that information can be obtained through the temporal velocity of the neighboring particles. This makes the method very flexible because we can actually apply external force and change the temporal velocity in the middle of the iteration process. When this happens, the pressure starts adapting to the new environment and finally converges to a different distribution from that which was previously assumed. This enables us to strongly couple the fluid simulation with different types of simulation, as long as these simulations only have velocity-based constraints.

We define velocity-based constraints here. A constraint in a system is referred to as velocity-based if it only requires the temporal velocity of the system to iteratively update its constraint force or impulse, and we apply the difference of it to update the temporal velocity. A velocity-based constraint must not require knowledge of the information on other constraints that cannot be obtained through the temporal velocity of the system. By this definition, we can say that each particle $i$ has a velocity-based constraint, the constraint force of which is $p_{i}^{t+1}$. We can update the constraint force and reapply this force to the temporal velocity without knowing the pressure of the other particles explicitly.

When the iteration process is terminated, we adopt the temporal velocity as the velocity in the next time step by ${\boldsymbol{u}}_{i}^{t+1}={\boldsymbol{u}}_{i}^{*}$ for each particle $i$. Then, we update the position of each particle $i$ by ${\boldsymbol{r}}_{i}^{t+1}={\boldsymbol{r}}_{i}^{t}+h{\boldsymbol{u}}_{i}^{t+1}$ to finish the time step.

2.6 Smoothing Pressure

As we will show in Section 5, the computed value of raw pressure $p$ contains a lot of noise, even though the fluid as a whole behaves very smoothly. However, we can actually obtain a smooth and sufficiently accurate pressure field by smoothing $p$ to some extent. A similar approach is presented by Kondo [37] to decrease the heavy noise of pressure. In [37], the virial theorem is used to make the smoothing process physically meaningful. Considering the balance of smoothness and locality, we define smoothed pressure $\tilde{p}_{i}$ for each particle $i$ by:

[TABLE]

which is equal to the weighted average of $p_{i}$ with weighting function $w_{\text{smooth}}$. We use the following weighting function for smoothing:

[TABLE]

which is proportional to a standard weighting function in the SPH method and has the same effective radius as the weighting function used in the main fluid computation.

3 Rigid-body Simulation

In this section, we introduce a rigid-body simulation method that can be strongly coupled with the incompressible fluid simulation introduced in the last section. As mentioned previously, in order to enable strong coupling, the rigid-body simulation must only have velocity-based constraints. This requirement, however, is not hard to satisfy because most of the constraints imposed on rigid bodies can be written as velocity-based constraints, which we defined in the last section.

We use ${\boldsymbol{x}}$, ${\boldsymbol{v}}$, ${\boldsymbol{\omega}}$, $m$, ${\boldsymbol{I}}$, ${\boldsymbol{R}}$, and ${\boldsymbol{q}}$ to represent a the center of gravity, linear velocity, angular velocity, mass, moment of inertia, rotation matrix, and a unit quaternion that represents the orientation, respectively, of a rigid body.

3.1 Time Integration

Let $h$ be the time step size. We use the same notation as in Section 2. We use $\phi^{t}$ to denote physical quantity $\phi$ at time $t$. We integrate the center of gravity ${\boldsymbol{x}}$ of a rigid body as follows:

[TABLE]

In addition, we integrate quaternion ${\boldsymbol{q}}$ of a rigid body as follows:

[TABLE]

where ${\boldsymbol{d}}{\boldsymbol{q}}^{t+1}$ is a unit quaternion that represents a rotation, the rotation axis of which is ${\boldsymbol{\omega}}^{t+1}$ and the rotation angle of which is $h{\left\lVert{\boldsymbol{\omega}}^{t+1}\right\rVert}$. Since the moment of inertia tensor ${\boldsymbol{I}}$ depends on the orientation of the rigid body, using only (45) does not accurately conserve angular momentum ${\boldsymbol{I}}{\boldsymbol{\omega}}$. However, computing an accurate angular velocity requires implicit calculation due to a numerical instability and is therefore costly. In the present paper, we assume that a change in angular velocity only occurs due to an external force applied to the rigid body. If the conservation of the angular momentum is important, a gyroscopic force can be calculated as an external force, which can be applied to a rigid body.

The linear velocity and the angular velocity at time $t+1$ is computed as follows. Starting from the previous velocities ${\boldsymbol{v}}^{t}$ and ${\boldsymbol{\omega}}^{t}$, we first compute temporal velocities ${\boldsymbol{v}}^{*}$ and ${\boldsymbol{\omega}}^{*}$ by applying external forces explicitly and then correct velocities by applying constraint impulses iteratively, so that they satisfy all of the velocity-based constraints. Finally, we adopt the corrected temporal velocities as the velocities at time $t+1$ as ${\boldsymbol{v}}^{t+1}={\boldsymbol{v}}^{*}$ and ${\boldsymbol{\omega}}^{t+1}={\boldsymbol{\omega}}^{*}$.

3.2 Contact Constraint

When two rigid bodies touch or collide, we generate contact points between these bodies and handle these bodies as contact constraints. Contact points are usually generated as they form the convex hull of the touching area. In the present paper, each contact point is treated as an independent constraint, and we compute the constraint impulse of each contact point independently.

Consider that two rigid bodies, namely, $A$ and $B$, are touching at point ${\boldsymbol{p}}$ with the unit contact normal ${\boldsymbol{n}}$, where the contact normal is directed from rigid body $B$ to rigid body $A$ at the contact point. We use subscript $\phi_{X}$ to denote arbitrary physical quantity $\phi$ of rigid body $X$. Let ${\boldsymbol{r}}_{A}={\boldsymbol{p}}-{\boldsymbol{x}}_{A}$ and ${\boldsymbol{r}}_{B}={\boldsymbol{p}}-{\boldsymbol{x}}_{B}$ be the relative contact positions, and let $N$ be the contact impulse along ${\boldsymbol{n}}$. Then, impulses $N{\boldsymbol{n}}$ and $-N{\boldsymbol{n}}$ are applied to rigid bodies $A$ and $B$, respectively. The changes in velocities of each rigid body caused by the impulses are as follows:

[TABLE]

We define the relative velocity at contact point ${\boldsymbol{v}}_{AB}$ by

[TABLE]

Then, the change in relative velocity $\Delta{\boldsymbol{v}}_{AB}$ can be written as

[TABLE]

where $[{\boldsymbol{a}}\times]$ is a skew-symmetric matrix that satisfies $[{\boldsymbol{a}}\times]{\boldsymbol{b}}={\boldsymbol{a}}\times{\boldsymbol{b}}$ for arbitrary vectors ${\boldsymbol{a}}$ and ${\boldsymbol{b}}$, and ${\boldsymbol{1}}$ is the identity matrix of order three. By taking the dot product of (56) and ${\boldsymbol{n}}$, we obtain the change in relative velocity along the normal:

[TABLE]

Here, ${\left\lVert{\boldsymbol{a}}\right\rVert}_{{\boldsymbol{A}}}=\sqrt{{\boldsymbol{a}}^{T}{\boldsymbol{A}}{\boldsymbol{a}}}$ for an arbitrary vector ${\boldsymbol{a}}$ and symmetric positive definite matrix ${\boldsymbol{A}}$, and the effective mass along normal $m_{AB{\boldsymbol{n}}}$ is defined as

[TABLE]

From (59), we obtain the linear relationship between $\Delta{\boldsymbol{v}}_{AB}\cdot{\boldsymbol{n}}$ and $N$. In order to change ${\boldsymbol{v}}_{AB}\cdot{\boldsymbol{n}}$ by some amount $\varepsilon$, we need to set $N$ to $m_{AB{\boldsymbol{n}}}\varepsilon$. Note that in actual computation, we iteratively update the value of $N$ as well as the temporal velocities of the rigid bodies.

The constraint of the contact point can be written as follows in the form of an LCP:

[TABLE]

where $b_{AB}$ is the target velocity of the contact point. The target velocity represents how fast the rigid bodies move away from each other along the normal after the collision. Using the coefficient of restitution $e$, we can estimate the relative velocity along the normal after the collision as $-e{\boldsymbol{v}}_{AB}^{*}\cdot{\boldsymbol{n}}$. Note that this value must be observed and fixed before any constraint impulse is applied. Although this estimation is not always physically correct because the energy passing through multiple bodies is not tracked, the estimation gives satisfactory results in most cases. Newton’s cradle is a typical case in which this estimation does not work well. Since the estimation depends only on the velocities, simply setting $b_{AB}=-e{\boldsymbol{v}}_{AB}^{*}\cdot{\boldsymbol{n}}$ will let rigid bodies gradually pass through each other due to numerical error if we do not correct positions separately. This is similar to the case of the incompressibility of fluid shown in Section 2.5. We should either correct positions or mix a positional term with the constraint in order to maintain the incompressibility. We also mix a positional term with the constraint in the case of the rigid-body simulation. Let $d_{AB}$ be the penetration depth of rigid bodies $A$ and $B$ at the contact point, and let $\alpha\in[0,1]$ be a constant to control the amount of mixing. Then, we set $b_{AB}$ as

[TABLE]

Contact constraints are velocity-based constraints, and we can solve these constraints in the same way as in Section 2.5. We iterate the following steps for each contact constraint, until some stopping criteria are met. First, we compute the difference $\delta_{AB}$ of the target velocity and the current relative velocity along the normal as

[TABLE]

Then, we compute the next value of impulse $N^{\prime}$ by

[TABLE]

Finally, we update impulse $N$ and apply its difference to each rigid body to update the temporal velocities:

[TABLE]

4 Rigid Body and Fluid Interaction

We introduced the incompressible fluid computation in Section 2, and the rigid-body computation in Section 3. In this section, we introduce the strongly coupled simulation method of rigid bodies and fluid.

We defined velocity-based constraints that can be strongly coupled with other velocity-based constraints. Thus, for the strongly coupled simulation, we should define the rigid-body and fluid interaction model, which consists of velocity-based constraints.

4.1 Shape Representation of a Rigid Body

In the present paper, we adopt particles to represent the shape of a rigid body. Since particles are used to represent fluid, there are many benefits to using particles as rigid bodies.

To construct a rigid body with particles, we first fill the inside of the rigid bodies with particles. Each particle is placed at the center of a grid of interval $l$, which is the same value as the interval of initial particle arrangement for fluid. Then, we compute the mass, the center of gravity, and the inertia tensor matrix of the rigid body. In this step, particles are considered to be axis-aligned boxes of edge length $l$. The position of each particle relative to the center of gravity of the rigid body is computed and stored in order to update the position of the particle when the position and orientation of the rigid body are changed.

During collision detection between rigid bodies, each particle is treated as a sphere of radius $l/2$. In addition, contact information is generated based on the sphere shape. In other words, when a pair of particles collides, a contact point is generated at the middle point between two particles and the contact normal is set to a vector parallel to the relative position vector of two particles. The penetration depth of the contact point is set to $l-r$, where $r$ is the distance between two particles.

In order to compute the particle number density of fluid particles around the boundary of a rigid body and fluid accurately, we include rigid-body particles in the particle number density computation of fluid particles. The particle number density of a fluid particle is computed using all kinds of particles.

4.2 Contact Constraint between a Rigid Body and a Fluid Particle

Since rigid bodies and fluid are both represented as particles, contact computation between rigid bodies and fluid can be performed in a particle-based manner. We treat fluid particles as rigid-body particles during the rigid-body and fluid interaction computation. Each fluid particle is considered as a rigid-body particle, which is an independent rigid body. The physical properties of the rigid body are defined based on a sphere of radius $l/2$ and mass $\rho l^{d}$, where $\rho$ is the density of the fluid and $d$ is the dimension of the simulation. The sphere is placed at the same position as the particle, and the rigid body shares its linear velocity with the velocity of the particle. That means that the linear velocity of the rigid body actually refers to the velocity of the particle, and if one velocity is changed, then the other velocity is also changed to the same value. This relationship is also applied to the temporal linear velocity of the rigid body and the temporal velocity of the particle while the iteration for the constraint resolution is running.

Collision detection and contact generation between rigid bodies and rigid bodies for fluid particles are performed in the same way as those described in the previous subsection. Since every interaction is treated as a rigid-body contact constraint, and a rigid-body contact constraint is a velocity-based constraint, we only need use velocity-based constraints for the entire simulation. This enables us to strongly couple rigid bodies and fluid in the simulation. We finally show the computation procedure of the entire time step for the strongly coupled simulation in Figure 1.

5 Computation Examples

In this section, we show computation examples of the proposed method. In Section 5.1, we compute a scene that consists of multiple rigid bodies in order to verify inter-rigid-body collision. In Section 5.2, Section 5.4, and Section 5.3, we compute scenes that consist of incompressible fluid in order to verify dynamic and static properties of incompressible fluid. In Section 5.5 and Section 5.6, we compute scenes that consist of both rigid bodies and incompressible fluid in order to verify interaction and strongly-coupled properties between rigid bodies and fluid. Finally, in Section 5.7, we provide a complex scene that involves multiple inter-rigid-body contacts and interaction between rigid bodies and fluid.

5.1 Rigid-Body Computation

We show a two-dimensional computation example of rigid-body interaction we introduced in Section 3. The initial configuration of the simulation is shown in Figure 2. Each box has the same mass, and the floor below the boxes has infinite mass and is not affected by gravity. The coefficient of restitution is set to 0.2 for all contact points. The interval of the particle arrangement $l$ is equal to 0.03 m, and the time step $h$ is 0.005 s. The gravitational acceleration $g$ is set to 9.8 m/s2 Snapshots of the simulation are shown in Figure 3. The piled boxes are dropped and collapsed, scattering over the floor. Thanks to the mixed term of constraint (64), no penetration between rigid bodies can be found in Figure 3. Although no friction is considered, we can see the friction-like effect in the simulation since the surfaces of the rigid bodies consist of particles and are uneven.

5.2 Hydrostatic Pressure Computation

Since we introduce a novel method for fluid computation, it is worthwhile to confirm the basic properties of the fluid. We compute the hydrostatic pressure of fluid. The initial configuration is shown in Figure 4. We set $l$ to 0.02 m, $r_{e}/l$ to 2.1, $g$ to 10.0 m/s2, and $h$ to 0.001 s. We show the distributions of the raw pressure $p$ and the smoothed pressure $\tilde{p}$ that is computed according to (42) at time 1.0 s in Figure 5. From Figure 5, we can observe that the raw pressure distribution suffers from heavy noise, whereas the smoothed distribution appears better and to be accurate. As we describe in Section 2.2, the smoothing length is equal to $r_{e}$, and thus the noise can be said to have only a local and limited effect. In addition, the distribution of the deviation from the standard particle number density $(n_{i}-n_{0})/n_{0}$ for each particle $i$ at time 1.0 s is shown in Figure 6. From the figure, we can see that no particle is compressed to the level of 1.0% deviation, and the incompressibility is highly assured.

5.3 Dam-break Computation

The dam-break problem is simulated using the proposed method. The initial configuration is shown in Figure 7. This is the same configuration as in the experiment performed by Koshizuka and Oka [23]. We set $l$ to 0.005 m, $r_{e}/l$ to 2.1, $g$ to 9.80665 m/s2, and $h$ to 0.0005 s. We show the snapshots of the simulation in Figure 8, which is colored according to smoothed pressure. We can observe that the fluid reaches the right wall approximately $t$ = 0.3 s after the beginning of the simulation and causes a heavy splash. The fluid makes a breaking wave from approximately $t$ = 0.7 s to $t$ = 0.9 s. The features of the simulation visibly match those of the experiment.

5.4 Circular Patch Computation

An extending circular droplet is computed using the proposed method. The radius of the droplet is set to 0.5 m, and the center of the droplet is located at the origin. The initial velocity of each particle $i$ at $(x,y)$ is given by

[TABLE]

We set $l$ to 0.005 m, $r_{e}/l$ to 2.5, $g$ to 0 m/s2, and $h$ to 0.002 s. This problem can be solved analytically [38] and the theoretical solution is given by

[TABLE]

with boundary conditions

[TABLE]

where $a$ and $b$ are the semi-minor axis and the semi-major axis of the droplet in meters, respectively. Snapshots are shown in Figure 9, and the simulation result is shown in Figure 10. The theoretical result in Figure 10 is computed numerically. Although pressure disturbance can be observed as the droplet extends due to the disorder of the particle arrangement, the graph shows the simulated values of the semi-minor axis and the semi-major axis well match the theoretical result.

5.5 Buoyancy Computation

In order to confirm that the interaction between a rigid body and fluid is computed correctly, we put boxes with different densities into a fluid to observe the effect of buoyancy. The initial configuration of the simulation is shown in Figure 11. The box is a square of edge length 0.6 m, and its density is set to have a specified density ratio to the fluid. The density ratio varies from 0.1 to 0.9 in steps of 0.1. We set $l$ to 0.01 m, $r_{e}/l$ to 2.1, $g$ to 9.80665 m/s2, and $h$ to 0.0025 s. We measure the ratio of the submerged volume of the box when the box becomes stationary in the fluid, and the result is shown in Figure 12. The graph shows that the simulated values match the theoretical values well.

5.6 Seesaw Computation

In strongly coupled simulations, multiple interactions of rigid bodies and fluid can be computed simultaneously. To confirm this, we run the following seesaw computation. The initial configuration of the simulation is shown in Figure 13. The densities of the rigid body and fluid are 490 kg/m2 and 1,000 kg/m2, respectively. The top square part of the fluid has an initial velocity of 1.0 m/s toward the left, and other fluid and the rigid body are set to be still. The rigid body is pinned at and can only rotate around its center of gravity. We set $l$ to 0.02 m, $r_{e}/l$ to 2.1, $g$ to 0.0 m/s2, and $h$ to 0.005 s. If the simulation is strongly coupled, all of the following things happen instantly in a single time step. The top part of the fluid pushes the rigid body toward the left so that the rigid body starts to rotate and pushes the bottom part of the fluid toward the right. This results in positive pressure at the bottom part of the fluid. Figure 14 shows the pressure distribution and the velocity distribution just after solving all constraints and before updating positions. We can observe that the rigid body is rotating and the bottom part of the fluid has positive pressure due to the rotation. We also show the result obtained with the PMS model in Figure 15. Since the PMS model cannot handle strong coupling, the bottom part of the fluid has no positive pressure and zero velocity. In subsequent steps, however, the PMS model can handle weakly-coupled interaction between the rigid body and the lower part of the fluid and positive pressure will be observed there. Note that the result obtained with the PMS model shows a lower peak pressure, as compared to the proposed method, because the top part of the fluid cannot take into account the existence of the bottom part of the fluid. As such, it is easier for the top part of the fluid to push the rigid body.

5.7 Complex Scene

Finally, we give a computation example of a complex and dynamic scene that contains both inter-rigid-body collision and rigid bodies and fluid interaction. The centers of gravity of two cross-shaped rigid bodies at the bottom are fixed and can only rotate around the centers of gravity. The other rigid bodies are not fixed and can move freely. The simulation result is shown in Figure 16. We set $l$ to 0.02 m, $r_{e}/l$ to 2.1, $g$ to 9.80665 m/s2, and $h$ to 0.001 s. We can observe that the simulation visibly runs without problems.

6 Conclusion

In this research, we proposed a method to simulate an incompressible fluid that uses an LCP to formulate the incompressibility constraint and enables strong coupling with rigid bodies. We formulated velocity-based constraints that generalize incompressibility constraints in fluid computation and non-penetration constraints in rigid-body computation, which provides the general framework of various strongly coupled simulations. Through numerical examples, we have demonstrated that the proposed method can compute incompressible fluid accurately and achieves strong coupling with rigid bodies correctly. With this method, we can use ordinary impulse-based methods for rigid-body simulation; therefore, the proposed method fits in well with existing software for rigid-body simulation.

The remaining problems of the proposed method include the fact that it is difficult to compute negative pressure because we required the pressure to be nonnegative in order not to cause an attractive force between particles. Some stabilization technique is likely needed to allow for negative pressure. In addition, generalization of the shape representation of rigid bodies is desirable so that we can represent a smoother face and decrease the computation time by means other than the use of particles. We leave these as future research topics.