# Rigid Body Sampling and Individual Time Stepping for Rigid-Fluid Coupling of Fluid Simulation

Authors: Xiaokun Wang; Xiaojuan Ban; Yalan Zhang; Xu Liu
Journal: Scientific Programming
 (2017)

---

## Abstract

In this paper, we propose an efficient and simple rigid-fluid coupling scheme with scientific programming algorithms for particle-based fluid simulation and three-dimensional visualization. Our approach samples the surface of rigid bodies with boundary particles that interact with fluids. It contains two procedures, that is, surface sampling and sampling relaxation, which insures uniform distribution of particles with less iterations. Furthermore, we present a rigid-fluid coupling scheme integrating individual time stepping to rigid-fluid coupling, which gains an obvious speedup compared to previous method. The experimental results demonstrate the effectiveness of our approach.


---

## Body

## 1. Introduction

Physically based fluid simulation is a popular issue in computer graphics and virtual reality while having a huge research and application demand in three-dimensional visualization and human- computer interactions. More realistic effects and higher simulation efficiency are the main goals; thus, reasonable, efficient, and scientific programming algorithms are needed to design and implement the animation. Two major schemes are employed for animating fluids: the grid-based Eulerian approach and particle-based Lagrangian approach. Eulerian method is particularly suited to simulate large volumes fluid, while being restricted by time step and computing time for small scale features. In contrast, Lagrangian method is suitable for capturing small scale effects such as spindrift and droplet. Among various particle-based approaches, Smoothed Particle Hydrodynamics (SPH) is the most popular method for simulating fluid due to computational simplicity and efficiency.In reality, rigid-fluid interaction widely exists in many scenarios. As a result, the interesting fluid behaviors emerge once rigid objects are added to fluid simulation. While interaction between particle-based fluids and rigid objects seems to be straightforward, there are still several issues not well resolved. For one thing, rigid bodies must be sampled to particles in order to interact with particle-based fluids, but only a few rigid boundary sampling methods can be directly employed in rigid-fluid coupling simulation. For another thing, the computational expenses of rigid-fluid coupling are considerable. To deal with the increasing demands for more detailed fluids and high efficiency, we present rigid sampling and individual time stepping for rigid-fluid coupling and design a practical and easy rigid-fluid animation simulation scheme with our scientific programming algorithms.


## 2. Related Work

Desbrun and Gascuel introduced SPH to computer graphics for simulating deformable objects [1]. SPH became popular in computer graphics for various fluid phenomena. Monaghan addressed simulating free surface flows with SPH [2] that serves as a basis for SPH fluid simulation. Muller et al. [3] proposed using gas state equation with surface tension and viscosity forces for fluid simulation, which also bring compressibility issue. Becker and Teschner [4] proposed WCSPH employing Tait equation to reduce compressibility. It significantly increased realistic effects but the efficiency is limited by time step. As incompressibility expenses computation time, many improved algorithms were addressed to enhance the efficiency. Solenthaler and Pajarola [5] presented PCISPH using a prediction-correction scheme to determine particle’s pressure and large time steps which significantly improved efficiency comparing to WCSPH. Another similar method that ensures incompressibility by iterative process is LPSPH [6]. Afterwards, Ihmsen et al. addressed a more efficient approach IISPH [7]. It carefully constructed pressure Poisson equation and solved the linear system using Relaxed Jacobi, which has a great improvement in stability as well as convergence speed and is particularly suitable for large-scale scene. Recently, a promising approach for impressible SPH has been proposed by Bender and Koschier [8]. It combines two pressure solvers which enforce low volume compression and a divergence-free velocity field and permits large time steps that yields a considerable performance gain.Besides, adaptive method, either spatial resolution or time discretization, is another way to promote efficiency. They allot computing resources to regions with complex flow behavior. Space adaptive methods [9–11] adaptively sample particle and employ less particles to produce similar details. Large particles are divided into small particles if high resolution is needed, and vice versa. However, difficulties exist in reproducing quantity when refining particles, and neighbor searching is usually the bottleneck. As an alternate to adaptive spatial discretization, the time domain can be adaptively sampled as well. Globally adaptive time stepping methods [1, 12, 13] employ a single time step to adjust each step for all particles in consideration of CFL condition. Though each particle has the current smallest time step, it is not the most efficient way. Locally adaptive time stepping methods [9, 14, 15] use different time steps for particles. Desbrun and Cani proposed that each particle evaluates forces depending on its current individual time step [9]. He et al. [16] adopt this idea and implement stable simulation of stiff fluids. It updates position, velocity, and density for active particles and interpolates for inactive particles. In this paper, we integrate it to rigid-fluid coupling to reduce the computational time.For boundary handling in SPH fluid simulation, distance-based penalty methods were commonly employed [17–19]. Nonetheless, these methods require large penalty forces which limit the time step and make particles stick to the boundary on account of lacking fluid neighbors. Frozen or ghost particles based models are used to solve the problem of sticking particles [20]. In order to avoid penetration, more than one layer of frozen particles were used [21], or the positions of penetrating particles should be corrected [13]. However, handling two-way interaction is troublesome since the elevated density near boundary in one phase affects particles in the other phase. For this reason and for the lack of fluid neighbors, Ghost SPH [22] solved those problem using a narrow layer of ghost particles and Akinci et al. employed boundary particles to correct the calculation of fluid density [23]. Because Ghost SPH is more time consuming, we use Akinci’s boundary handling method which is simple and easy to achieve in this paper.For rigid-fluid coupling, several approaches were presented up to now. In [24], the fluid is represented as rigid spheres and switching impulses with rigid bodies. In [25, 26], the pressure at the boundary is taken into consideration for two-way fluid-rigid coupling. Then, Oh et al. proposed an impulse-based scheme for two-way coupling of SPH fluids with rigid bodies [27]. Becker et al. presented direct forcing for rigid-fluid coupling [28] which employs a prediction-correction scheme to enforce particle positions and velocities to specific values. Akinci et al. presented a momentum-conserving two-way coupling approach based on hydrodynamic forces that use boundary particles to sample the surface of rigid bodies [23]. We present a rigid-fluid coupling scheme by integrating individual time stepping to Akinci’s boundary handling method that gains an obvious speedup [29].Rigid bodies sampling includes particle-based methods and polygonization-based methods. Turk repelled particles on surfaces to get a uniform sample [30] and also simplified a polygonization through reducing the number of polygons [31]. Witkin and Heckbert employed local repulsion to make particles spread uniform [32]. Nehab and Shilane [33] presented algorithm of stratified point sampling. Cook addressed stochastic sampling of Poisson disk distributions with blue noise [34]. Blue sampling has the ability to generate random points and get uniform distribution of sampling points set. Therefore, the following sampling methods always have blue noise characteristics. Corsini et al. sampled triangular meshes with blue noise properties [35]. Dunbar and Humphreys [36] modified Poisson disk sample using a spatial data structure. Bridson [37] simplified Dunbar’s approach with rejection sampling and extending it to higher dimensions. Then, Schechter [22] modified Bridson’s approach and employed it to Ghost SPH. Inspired by Schechter’s approach, we address sampling method improved by SPH equation that is more efficient and easy to implement.


## 3. Particle-Based Fluid Simulation Framework

In particle-based fluid simulation, the forces acting on particles are derived from the Navier-Stokes equations. The conservation of mass and momentum are written as(1)dρidt=-ρi∇·vi,ρiDviDt=-∇pi+ρig+μ∇2vi,where vi is the velocity, ρi is the density, pi is the pressure, μ is the viscosity coefficient, and g is the external force field.SPH works by obtaining approximate numerical solutions of fluid dynamics equations by expressing fluids with particles. In SPH, the representation of a field variableA at location xi is defined as(2)Axi=∑jmjAjρjWxi-xj,h,where mj and ρj represent particle mass and density, respectively, Wxi-xj,h is smoothing kernel, and h is smoothing radius.It can be easily derived from the basic SPH equation by substituting fluid densityρ into (2), that is,(3)ρi=∑jmjWxi-xj,h.Therefore, particles’ pressure force fiP and viscous force fiv can be written as(4)fiP=-∑jmjPiρi2+Pjρj2∇Wij,fiv=μ∑jmjvjiρj∇2Wij.In this paper, we use Tait equation [4] to calculate the pressure; that is, pi=(ρ0cS2/γ)ρi/ρ0γ-1, where ρ0=1000 is the rest density of the fluid, γ=7 is stiffness parameter, and cS is velocity of sound. We use the equation in [23] to compute viscous force.


## 4. Boundary Handling for Particle-Based Fluids

In rigid-fluid interaction, there are three types of forces among particles: the forces between fluid particles and boundary particlesForcerigid-fluid, the forces between fluid particles Forcefluid-fluid, the forces between boundary particles Forcerigid-rigid, which is shown in Figure 1. In our simulation, we sample rigid body to obtain boundary particles which is described in Section 5. Moreover, for the boundary handling way of rigid-fluid interaction, we implement our simulation based on the work of [23]. The following will briefly introduce boundary handling in this section.Figure 1
Forces between boundary particles and fluid particles.Considering influence of boundary particles, density formula of fluid particles in (3) needs to introduce weighted summation influence of boundary particle [23], that is,(5)ρfi=∑jmfjWij+∑kmbkWik,where fj and bk denote fluid particle j and boundary particle k, respectively. The first summation calculates the affection of adjacent fluid particles, while the second summation computes the influence of adjacent boundary particles. This formula can overcome the problem of boundary defects in SPH fluid simulation to some extent.Due to the use of boundary particle mass in (5), the density of fluid particles is incorrect or unstable when the boundary particle density is set unreasonably or unevenly distributed. Hence, consider the contribution of boundary particles to fluid particle by the volume of boundary particles is(6)Ψbiρ0=ρ0Vbiwhere ρ0 denotes the remaining density of fluid and Vbi is the estimation value of boundary area volume of corresponding boundary particles. Applying Ψbiρ0 to replace the boundary particle mass can guarantee the stability.Thus, (5) can be rewritten as (7)ρfi=∑jmfjWij+∑kΨbkρ0iWik.The most important interaction between fluid particles and boundary particles is the pressure. The pressure acceleration generated by boundary particles to fluid particles can be computed as(8)dvfidt=-kpfiρfi2∑kΨbkρ0i∇Wik,where pfi>0 takes k=2. When pfi<0, boundary particles and fluid particles attract each other; then, we can adjust parameter k  (0≤k≤2) to realize different adsorption effects, and we choose k=1 in our experiment.To simulate the friction between fluid and container wall or the interaction of rigid body and fluid, we have to compute the friction of boundary particles with fluid particles. The friction is calculated by the artificial viscosity; that is,(9)dvfidt=-∑kΨbkρ0iΠik∇Wik,where Πik=-νvikTxik/(xik2+εh2), ν=2αhcs/(ρk+ρj).On account of (8) and (9) that listed the forces for fluid particles, we can get the forces of boundary particles using Newton’s third law. The forces generated by fluid particles to boundary particles are(10)Fbk=∑ikpfiρfi2+ΠikmfiΨbkρ0i∇Wik,where i denotes the fluid neighbors of boundary particle k. It is the counteracting force of (8) and (9).For a rigid body, the total force and torque need to be calculated. This can be separately written as(11)Frigid=∑kFbk,τrigid=∑kxk-xrigidcm×Fbk,where xk is the location of boundary particle k and xrigidcm denotes the mass center of a rigid body. The total force and torque will be transmitted to the physics engine to handle the motion of rigid bodies.


## 5. Rigid Boundary Sampling

Rigid body sampling is the first issue in rigid-fluid coupling which we have to handle. We propose a rigid body sampling algorithm which is an extension of Poisson disk method and sampling method in [22] for rigid-fluid coupling.For rigid objects sampling, boundary particles are used to sample the surface of rigid bodies that has several merits. For one thing, using particles permits us to get a rigid model that can handle different shapes even with complex geometry structure. For another thing, the use of boundary particles successfully alleviates sticking artifacts and makes sampling uniform.There are two components in our sampling: surface sampling and surface relaxation. As shown in Figure2, it first samples the surface of rigid object image, and then it improves initial sampling with surface relaxation. In order to realize the first procedure, it needs fast projection of points to the surface. Hence, level set method is employed to express surface geometry with ϕ>0 and ϕ<0 denoting exterior and interior of rigid objects, respectively, while ϕ=0 denotes surface of rigid objects.Figure 2
Surface sampling and relaxation. (a) Surface sampling. (b) Surface relaxation. Black points: newly added sample. Gray points: surface sampling points. White points: exterior points before being projected to the surface.
(a)
(b)After obtaining the surface geometry of signed distance function, we use surface sampling method proposed in [22] (as shown in Algorithm 1), first, searching seed points on the surface by checking each grid cell intersecting with the surface; the details are as follows: projecting random points from the cell to the surface and stopping when the point satisfies the Poisson disk criterion, which operates k attempts in a cell; when obtaining a seed sample, continuing to sample it and taking a step of size e·r from the previous samples along a random tangential direction d; then projecting to the surface and checking the Poisson disk criterion again. Parameters were chosen as k=30 and e=1.085.Algorithm 1:Surface sampling.
Input: Level set ϕ, radius r, count t, constant e
Output:Sample set S
(1)forall grid cells C that ϕ changes signdo
(2)foreach tdo
(3)    Generate random pointpin C
(4)    Projectpto surface of ϕ
(5)ifpsatisfies Poisson Disk criterion in Sthen
(6)S←S∪{p}
(7)Break
(8)ifno point was found in Cthen
(9)Continue
(10)whilenew samples are founddo
(11)    Generate random tangential directiondof surface atp
(12)q←p+d⋅e⋅r
(13)Projectqto surface of ϕ
(14)    ifq satisfies the Poisson Disk criterion in Sthen
(15)S←S∪{q}
(16)p←qIn order to optimize the position of sample points, reduce noises, and get a uniform distribution set of sampling points, we need to further improve sampling set with a surface relaxation step. Inspired by the relaxation algorithm proposed in [22] and SPH interpolation method, surface relaxation algorithm is presented in Algorithm 2. Unlike employing random testing way in [22], we compel particles to move according to the density gradient. This ensures that the particles move to sparse place, so as to insure uniform distribution of particles.Algorithm 2:Modified surface relaxation.
Input: sample setS, Level setϕ, radiusr, count t, constant f
Output:relaxed sample set S
(1)foreachtdo
(2)foreachpi∈Sdo
(3)     compute densityρi(t), average densityρ-(t)
(4)     compute density gradient∇ρi(t)
(5)d←r⋅(ρit-ρ-t)/ρ-(t)⋅f
(6)pnew←p+d⋅∇ρi(t)
(7)      ifpnewoutsideϕor came from surface sample
(8)Projectpnewto surface of ϕ
(9)     ifpnewsatisfies the Poisson Disk criterion in S
(10)p←pnewIt starts with the initial particles seed in Algorithm1 and attempts to reposition each sample through density gradient. Next it computes density ρi(t) and density gradient ∇ρi(t) of each surface particles and employs deviation of density ρi(t) and average density ρi(t)¯ as a coefficient to tune distance d. Then, it employs d·∇ρi(t) to adjust particle locations. Surface sample candidates are projected to the surface once again and are reserved which satisfies the Poisson disk criterion. Parametert is iterations andf is distance coefficient. According to SPH gradient formula, particle’s density gradient can be written as(12)∇ρi=∑j=1Nmj∇Wxi-xj,h.On the basis of signed distance field ϕ, it is fairly convenient for calculation of steps (7) and (8). If ϕ(pnew)≠0, it means pnew is not on the surface. While for projecting particles to surface, we compute signed distance field gradient ∇ϕ(pnew). Projection formula is (13)pnew=p+d·∇ρit-∇ϕpnew.We have done the experiment in 2 dimensions contrast to the relaxation algorithm in [22]. We randomly generate 100 points in a 0.1×0.1 square shown in Figure 3. The red points represent that it does not meet the conditions of Poisson disk.Figure 3
Initial samples in a0.1×0.1 square.Figure4 shows the relaxation results of Figure 3; the first row is our method and the second row is the method in [22]. The column (a) reveals the distribution of points after relaxation of two algorithms and the red points mean it does not satisfy Poisson disk condition. Each algorithm iterates 100 times, respectively, while column (b) illustrates the number of points that do not satisfy Poisson disk conditions for each iteration. It is obvious that our method can get a better effect with a slight concussion. Compared to the method in [22], the optimization effect is basically the same after 30 times’ iteration using our method while it is the optimal results of fast Poisson disk method. In addition, our method is more efficient. In MATLAB environment, all parameters are the same as mentioned in [22]; our method takes 2.44691 s while relaxation method in [22] costs 56.44153 s for 100 iterations.Figure 4
Relaxation results comparison of our method with fast Poisson disk. (a) Relaxation results. (b) The relationship between iteration times and points not conforming to the conditions.


## 6. Individual Time Stepping for Rigid-Fluid Coupling

In this section, we propose a rigid-fluid coupling method employing individual time stepping. As a result, larger time steps can be used comparing to previous methods and the overall computation time is reduced. In particle-based fluids, particles only interact with their neighbors. So permitting particles to have different time steps is more efficient than using a global time step for all particles. The individual time stepping computes time step for each particle and updates time step asynchronously.


### 6.1. Time Steps

For particle-based fluids, the time step must satisfies Courant-Friedrich-Levy (CFL) condition for numerical stability, that is,(14)ΔtCFL≤λvhvmax,where vmax=maxi⁡vi is the maximum velocity of particles and coefficient λv<1. In addition, it also has to consider particles’ maximum acceleration. Thus, the time step must also meet the condition(15)Δtf≤λfhfmax,where fmax=maxi⁡dvi/dt denotes the maximum force per unit mass of particles and λf<1. In [13], λv=0.4 and λf=0.25 are used for PCISPH, while we use λv=0.1 and λf=0.05 for WCSPH. Instead of using a constant time step, we adjust time step dynamically as(16)Δt≤minΔtCFL,Δtf.Thus, the time step for each particle i is(17)Δti=minjλvhvj,λfhdvj/dt,where j denotes that it iterates all the neighbors.However, the presented algorithm demands small coefficients for asynchrony, so we chooseλv=0.05 and λf=0.025.


### 6.2. Asynchronism

Since each particle has individual time step, we enforce asynchronous time integration for the update. To save computing resources, the system time step is chosen as the minimized individual time step:(18)Δt=mini⁡Δti,where Δti is computed by (17).The particlei will be updated if it satisfies the condition(19)tilast+Δti<t,where tilast denotes the last update time of particle i and t is the system time. It means that if system time is larger than the individual time step, particle i will be set as active particle and be updated.Semi-implicit Euler integration is generally used in SPH simulation. To accommodate for asynchronism, the semi-implicit Euler integrations can be expressed as(20)vitilast+Δt′=vitilast+Δt′aitilast,xitilast+Δt′=xitilast+Δt′vitilast+Δt′,where Δt′ is an independent integral time step which is different from the global time step Δt. For inactive particles, (20) is equivalent to interpolation, while, for active particles, they are semi-implicit Euler integration equations.


### 6.3. Algorithm

The individual time stepping for rigid-fluid coupling algorithm is shown in Algorithm3. In this algorithm, particle i has several extra variables; that is, dρit/dt denotes density derivative, Δti is time step, Δti′ is individual condition time step, and tilast is last updated time. In the algorithm, t is the system time and Δt is system update time step. Particle i is active if tilast+Δti<t.Algorithm 3:Individual time stepping for rigid-fluid coupling.
(1)whileanimatingdo
(2)   select active
(3)foreachactive fluid particle ido
(4)     find fluid and boundary neighbors
(5)for each fluid particle ido
(6)ifactivethen
(7)      computeρit, pit
(8)else
(9)      interpolateρit, pit using dρitilast/dt
(10)foreachactive fluid particle ido
(11)compute dvit/dt
(12)compute dρit/dt
(13)compute time step condition Δti′ (Eq. (15))
(14) tilast=t
(15)foreachboundary particle kdo
(16)compute forces (Eq. (10))
(17)foreachfluid particle ido
(18)    compute time stepΔti=minj⁡Δti′
(19)Δt=mini⁡(Δti)
(20)foreachrigid bodydo
(21)compute total forces, torques (Eq. (11))
(22)pass forces and torques to physics engine
(23)update rigid body
(24)update boundary particles of rigid body
(25)foreachfluid particle ido
(26) Δt′=t+Δt-tilast
(27) vitilast+Δt′=vitilast+Δt′(dvitilast/dt)
(28) xitilast+Δt′=xitilast+Δt′vitilast+Δt′
(29)t=t+ΔtIn order to analyse the proposed algorithm, we implement the breaking dam with obstacles experiment. The setting of this experiment is shown in Table1.Table 1
The setting of breaking dam with obstacles.
Item
Value
Simulation domain size
12
m
×
12
m
×
8
m
Fluid particles
153600
Boundary particles
73585
Smoothing kernel function
Cubic splines
Smoothing radius
0.2 m
Fluid particle width
0.1 mWe compare individual stepping method to adaptive stepping and constant stepping method in breaking dam with obstacles scene. The rending results are shown in Figure5 and the time statistics are listed in Table 2. From Figure 5, the fluid simulation results are almost not different using three methods, while in Table 2, we can find that our method gains 1.5 and 6.4 times speedup comparing to globally adaptive stepping method and constant stepping method, respectively. In addition, the average active particles percent of individual stepping method is 31%.Table 2
Comparison of experimental results of breaking dam with obstacles.
Method
Total comp. time
Avg.Δt (avg. active pct)
Speedup
Constant steps
175 min
0.11 ms
—
Globally adaptive
41 min
0.46 ms
—
Individual stepping
27 min
0.23 ms (31%)
1.5 (6.4)Figure 5
Rendering results of breaking dam with obstacles. (a) Individual stepping method. (b) Globally adaptive stepping method. (c) Constant stepping method.
(a)
(b)
(c)


## 6.1. Time Steps

For particle-based fluids, the time step must satisfies Courant-Friedrich-Levy (CFL) condition for numerical stability, that is,(14)ΔtCFL≤λvhvmax,where vmax=maxi⁡vi is the maximum velocity of particles and coefficient λv<1. In addition, it also has to consider particles’ maximum acceleration. Thus, the time step must also meet the condition(15)Δtf≤λfhfmax,where fmax=maxi⁡dvi/dt denotes the maximum force per unit mass of particles and λf<1. In [13], λv=0.4 and λf=0.25 are used for PCISPH, while we use λv=0.1 and λf=0.05 for WCSPH. Instead of using a constant time step, we adjust time step dynamically as(16)Δt≤minΔtCFL,Δtf.Thus, the time step for each particle i is(17)Δti=minjλvhvj,λfhdvj/dt,where j denotes that it iterates all the neighbors.However, the presented algorithm demands small coefficients for asynchrony, so we chooseλv=0.05 and λf=0.025.


## 6.2. Asynchronism

Since each particle has individual time step, we enforce asynchronous time integration for the update. To save computing resources, the system time step is chosen as the minimized individual time step:(18)Δt=mini⁡Δti,where Δti is computed by (17).The particlei will be updated if it satisfies the condition(19)tilast+Δti<t,where tilast denotes the last update time of particle i and t is the system time. It means that if system time is larger than the individual time step, particle i will be set as active particle and be updated.Semi-implicit Euler integration is generally used in SPH simulation. To accommodate for asynchronism, the semi-implicit Euler integrations can be expressed as(20)vitilast+Δt′=vitilast+Δt′aitilast,xitilast+Δt′=xitilast+Δt′vitilast+Δt′,where Δt′ is an independent integral time step which is different from the global time step Δt. For inactive particles, (20) is equivalent to interpolation, while, for active particles, they are semi-implicit Euler integration equations.


## 6.3. Algorithm

The individual time stepping for rigid-fluid coupling algorithm is shown in Algorithm3. In this algorithm, particle i has several extra variables; that is, dρit/dt denotes density derivative, Δti is time step, Δti′ is individual condition time step, and tilast is last updated time. In the algorithm, t is the system time and Δt is system update time step. Particle i is active if tilast+Δti<t.Algorithm 3:Individual time stepping for rigid-fluid coupling.
(1)whileanimatingdo
(2)   select active
(3)foreachactive fluid particle ido
(4)     find fluid and boundary neighbors
(5)for each fluid particle ido
(6)ifactivethen
(7)      computeρit, pit
(8)else
(9)      interpolateρit, pit using dρitilast/dt
(10)foreachactive fluid particle ido
(11)compute dvit/dt
(12)compute dρit/dt
(13)compute time step condition Δti′ (Eq. (15))
(14) tilast=t
(15)foreachboundary particle kdo
(16)compute forces (Eq. (10))
(17)foreachfluid particle ido
(18)    compute time stepΔti=minj⁡Δti′
(19)Δt=mini⁡(Δti)
(20)foreachrigid bodydo
(21)compute total forces, torques (Eq. (11))
(22)pass forces and torques to physics engine
(23)update rigid body
(24)update boundary particles of rigid body
(25)foreachfluid particle ido
(26) Δt′=t+Δt-tilast
(27) vitilast+Δt′=vitilast+Δt′(dvitilast/dt)
(28) xitilast+Δt′=xitilast+Δt′vitilast+Δt′
(29)t=t+ΔtIn order to analyse the proposed algorithm, we implement the breaking dam with obstacles experiment. The setting of this experiment is shown in Table1.Table 1
The setting of breaking dam with obstacles.
Item
Value
Simulation domain size
12
m
×
12
m
×
8
m
Fluid particles
153600
Boundary particles
73585
Smoothing kernel function
Cubic splines
Smoothing radius
0.2 m
Fluid particle width
0.1 mWe compare individual stepping method to adaptive stepping and constant stepping method in breaking dam with obstacles scene. The rending results are shown in Figure5 and the time statistics are listed in Table 2. From Figure 5, the fluid simulation results are almost not different using three methods, while in Table 2, we can find that our method gains 1.5 and 6.4 times speedup comparing to globally adaptive stepping method and constant stepping method, respectively. In addition, the average active particles percent of individual stepping method is 31%.Table 2
Comparison of experimental results of breaking dam with obstacles.
Method
Total comp. time
Avg.Δt (avg. active pct)
Speedup
Constant steps
175 min
0.11 ms
—
Globally adaptive
41 min
0.46 ms
—
Individual stepping
27 min
0.23 ms (31%)
1.5 (6.4)Figure 5
Rendering results of breaking dam with obstacles. (a) Individual stepping method. (b) Globally adaptive stepping method. (c) Constant stepping method.
(a)
(b)
(c)


## 7. Implementation and Results

All the experiments in this paper are implemented on Intel 3.50 GHz CPU with 4 cores. The simulation algorithms (Algorithms1, 2, and 3) and surface reconstruction [38, 39] are actualized with C++ language and multithreading technology. Bullet is used to simulate rigid objects while OpenMP served as parallelization. Images were rendered with Blender.To implement fluid-rigid coupling animation efficiently and scientifically, we design a fluid-rigid coupling programming simulation scheme shown in Figure6. We firstly initialize rigid and fluid particles configuration. Then, we do neighbor search using spatial hashing algorithm for each particle. Next, the governing equations of fluid, rigid boundary sampling, and boundary handling are solved as described in the previous sections. Then, total force and torque of rigid are calculated and provided to Bullet. After total force acting on particles is computed, particles are integrated to next time step using asynchronous update scheme introduced in Section 6.Figure 6
Flowchart of the programming scheme.In order to demonstrate the validity of the entire fluid-rigid coupling simulation system, we designed a scene of dropping multiple small squares into water. The setting and statistics are shown in Table3, and the experimental results are displayed in Figure 7. It can be seen from the diagram that small squares fall into water and splash water while they are rotated and inclined by water. Finally, small squares are force balance and floating on the water.Table 3
The setting and statistics of fluid-rigid coupling.
Item
Cubes fall into water
Water shock sculpture
Simulation domain size
12
m
×
12
m
×
12
m
10
m
×
8
m
×
5
m
Fluid particles
320 k
873 K
Boundary particles
47 K
156 K
Smoothing kernel function
Cubic splines
Cubic splines
Smoothing radius
0.05
0.1
Artificial viscosity coefficient
0.05
0.05
Surface tension coefficient
0
0
Rigid body mass
65 kg
1100 kg
Rigid body volume
0.125m-3
18.4 (4×2×2.3) m-3
Fluid rest density
1000 kg⋅m−3
1000 kg⋅m−3
SPH computing time (1 frame)
0.9008 min
8.46 min
Surface reconstruction time (1 frame)
74.0 s
261 s
Rendering time (1 frame)
13.85 min
8.86 minFigure 7
16 cubes fall into water. (a) Simulation in particle view; (b) rendering results.
(a)
(b)We realize another fluid-rigid coupling experiment displayed in Figure8. Figure 8(a) is the results in particle view while Figure 8(b) is the rendering results. The setting and statistics are illustrated in Table 3. In this scenario, the breaking dam of water hits the sculpture which is knocked down and pushed for some distances due to kinetic energy of water. The motions of sculpture are in line with expectations which proved that the simulation and calculation of fluid-rigid coupling system accord with physics laws.Figure 8
Water shock sculpture. (a) Simulation in particle view; (b) rendering results.
(a)
(b)This experiment further proved that our method can implement vivid fluid-rigid coupling animation simulation system with high realistic effects. It can be expected that this animation system can be used in virtual reality domain and special effects in film and game.


## 8. Conclusion

We proposed an efficient and simple rigid-fluid coupling scheme for particle-based fluid simulation. It samples rigid surface with boundary particles which are used to interact with fluids. It insures uniform distribution of particles which requires less iterations. In addition, we present an efficient rigid-fluid coupling approach combining individual time stepping with rigid-fluid coupling. Neighbors and forces of particles are updated only when needed, while computing resources are allocated to complex regions. It obtains an obvious speedup compared to previous methods. Besides, this scheme was integrated with rigid body coupling simulation with several scenes which has a good sense of visual reality. Overall, our method is efficient to compute while the sampling and coupling algorithm can be applied to other particle-based simulation or relevant approaches. Future work would be extending the proposed method to IISPH [7] or DFSPH [8] as well as large-scale scenarios.



---

*Source: 8502691-2017-03-02.xml*