# Real-Time Simulation of Fluid Scenes by Smoothed Particle Hydrodynamics and Marching Cubes

Authors: Weihong Wang; Zhongzhou Jiang; Honglin Qiu; Wei Li
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Simulating fluid scenes in 3DGIS is of great value in both theoretical research and practical applications. To achieve this goal, we present an algorithm for simulation of fluid scenes based on smoothed particle hydrodynamics. A 3D spatial grid partition algorithm is proposed to increase the speed for searching neighboring particles. We also propose a real-time interactive algorithm about particle and surface topography. We use Marching Cubes algorithm to extract the surface of free moving fluids from particles data. Experiments show that the algorithms improve the rate of rendering frame in realtime, reduce the computing time, and extract good real effects of fluid surface.


---

## Body

## 1. Introduction

Simulation of fluid scenes has a wide range of applications in the movie special effects, computer animation, game production, virtual reality, and many other fields [1–4]. At first, most studies of fluid simulation are based on the nonphysical procedural modeling method which generates scene by parameterized surface. Although the simulation is faster, these methods are not based on physical principles, thus lack of authenticity, and can’t simulate some detailed effects, such as the effect of wave overturning. At the same time these methods are based on random functions, it is difficult to achieve solid-liquid interaction. So most researchers focused on fluid simulation based on physical methods.At present, physics-based fluid simulation methods are mainly divided into two types [5]: one is grid-based Euler method, the other is particle-based Lagrange method. Smoothed particle hydrodynamics (SPH) is a new method based on particle-based Lagrange method, which is mainly used in astrophysics at first. J. Monaghan is the first one to apply SPH into simulation of free surface flow. Stam and Fiume [6] brought in SPH method into fluid simulation to achieve the effect of gas and flame. And then, the researches on fluid simulation mainly adopt particle-based Lagrange method. Harada et al. [7] made the searching adjacent particle feasible by constructing uniform spatial grid. Under this circumstances, fluid particles are relatively dispersed and lots of idle grid units would appear. Therefore, Grahn [8] made the fast searching adjacent particle of arbitrary size scene possible by constructing the space uniform grids with Hash function through GPU. Kolb and Cuntz [9] made use of GPU to achieve the whole process of SPH. However, this method calculates in the grids and interpolates in the particles, which leads to physical discontinuousness and simulation results distorted. David Lopez et al. [10] made use of SPH to build the physical pressure modal on river basin of the Villar Del Rey Dam, and it simulated the fluid scene of discharge floodwater. The way to deal with boundary interaction is too simple. So it simplified the terrain boundary to regular area, but the terrain surface is rugged in 3DGIS.Although there are many SPH methods to simulate fluid in literature, they can’t apply to simulate fluid in real-time and high efficiency in complex scene [11–13], such as simulating debris flow and flood and other process of fluid evolution in the 3DGIS. This paper improves 3D spatial grid partition algorithm to increase speed of neighboring particles searching, and we also propose a real-time interactive algorithm on particle and terrain surface. The results show that this method can calculate in real-time and has a good rendering rate.


## 2. Smoothed Particle Hydrodynamics

The basic idea of SPH makes fluid as a series of discrete particles, through the role of the smooth, which has certain radius of kernel functions to a particle physics scalar (such as density, pressure, etc.) assigning to the adjacent particles, as shown in Figure1. Physical scalar A of any particle X(x,y,z) can be calculated by
(2.1)A(x)=∑jmjAjρjW(X-Xj,h),(2.2)∇A(x)=∑jmjAjρj∇W(X-Xj,h),
where ρj is density of the particle Xj, mj is quality of the particle Xj, and h is radius of smoothing kernel. If the distance from the particle X to particle Xj, that is, |X-Xj|<h, the particles Xj is in smooth domain of the particle X, we can get the weight of particle Xj by  smoothing kernel function W(X). Otherwise |X-Xj|<h, the weight of the particles Xj is zero. In the solution of fluid equations, we often need to take derivatives of the physical scalar, and this operation only affects smoothing kernel function in SPH method. Therefore gradient of physical scalar value A can be expressed as (2.2).Figure 1
The basic principle of SPH.


## 3. Grid-Based Neighboring Particles Search

The existing neighboring particle search algorithms mostly adopt to traverse all particles directly, with time complexity ofO(n2) and n stands for the number of fluid particles. Chen et al. [14] constructed an index table of two-dimensional array of particle and spatial grid by dividing the three-dimensional space grid with grid number being the primary key and the index table being radix sorted (algorithm complexity of radix sort is O(n)). The index numbers of the first particle which corresponds to each grid are obtained so that each grid in the particle is identified.In this paper, the improved neighboring particle search algorithm doesn’t need to create the index table of particle and grid but divide grids, which saves memory space. It also doesn’t need to sort the particles and the index table, so it reduces the complexity of algorithm. The algorithm steps are as follows.Step 1.
Divide the three-dimensional spatial grid. In this paper, bounding box of the three-dimensional terrain space is divided into cube grid with length ofh. h stands for the radius of smoothing kernel. The Standard Template Library std::vector <Particle>m_Grid is used to store information of particles of each spatial grid. Particle is a class which stores density, velocity, spatial coordinates, and other physical information of fluid particle, and a Particle* pointer pointing to next particle.Step 2.
Put the fluid particles into three-dimensional spatial grid. Insert all the particlesXi(x,y,z)(1≤i≤n) into the corresponding spatial grid number (xi,yi,zi)  by (3.1). That is to insert the current particle pointer into the spatial grid:
(3.1)Grid(xi,yi,zi)={min_int(x/h),min_int(y/h),min_int(z/h)},
where x,y,z stands for the particle Xi in x,y,z axis coordinate, respectively, min_int is the positive integer function that get the minimum value, and grid(xi,yi,zi) represents spatial grid number.Step 3.
Calculate the neighboring particles of each particle. In SPH method, particles are only affected by other particles in radius of smoothing kernel, so we will find the scope for particle located in the neighborhood of 27 spatial grids. Firstly, calculate the spatial grid numberGrid(xi,yi,zi) of particle Xi by (3.1); then calculate index number of Grid(xi,yi,zi) by
(3.2)Gindex=yi*GridWidth*GridHeight+zi*GridWidth+xi.
GridWidth, GridDepth, and GridHeight, respectively, stand for the dimension on direction X,Y,Z of spatial grid, then the contained particles of m_Grid(Gindex) are stored in pointer array. Finally, visit all the particles in Grid{xi±1,yi±1,zi±1} and calculate current particle of density and physical pressure and so on by using the equations derivated from (2.1) and (2.2).


## 4. Interaction between the Fluid Particles and the Boundary

The simplest algorithm is exhaustive algorithm on dealing with interaction between fluids and boundary. Assuming that there areM geometric facets and N fluid particles in the scene, its time complexity is O(M*N). In the 3DGIS, the number of triangular patches of terrain mostly ranges from 105 to 107 and the number of fluid particles is generally between 103 and 106. If the exhaustive algorithm is adopted, the computation will be too huge to be accepted. In this regard, [14] proposed a boundary interaction algorithm based on spatial mesh, which inserts information of the terrain boundary and obstacle into the spatial grid, then judges whether the fluid particles are intersected with boundary geometry patches and obstacle in spatial grid when the fluid particles pass through. The algorithm requires geometry patches of boundaries to be small enough; if geometry patches are large, the number of grid number each grid space occupied will increase, which will lead to reduced efficiency of algorithm.Because the number of triangular grids of terrain is huge in the 3DGIS, if the above algorithms are adopted for boundary interaction, the amount of calculation will be too huge and the efficiency of real-time rendering will be degenerated. Regarding this, this paper proposes a real-time interaction method of the boundary bounce particles, where the time complexity of algorithms isO(N) and N is the number of fluid particles. Interaction between the fluid particles and the terrain is mainly used in interaction between the fluid particles and the terrain surface to prevent fluid particles from penetrating through the terrain surface, shown in Figure 2 The basic idea of algorithm is as follows.Figure 2
Interaction between the fluid particles and the terrain surface.Step 1.
Calculate intersectionTi of the ray and triangular grid of the terrain surface. The origin of the rays is Xi(xi,yi,zi) and direction is straight down (-y). As DEM data is regular grid data, two-dimension array can be used for storage. Spatial coordinate Vex1 of the triangle patches of terrain surface can be calculated by (4.1). Ver2 is (m*x_step,m_height[m][n+1],(n+1)*x_step). In a similar way, coordinates of Ver3 and Ver4 can be calculated. The x and z axis’s coordinates of Ti and Xi are equal, if (zi-z)/(xi-x)≥1, then calculate the elevation of Ti by bilinear interpolation of triangle 1. Otherwise, calculate the elevation of Ti by bilinear interpolation of triangle 2, xi and zi are the coordinates of Vi, x and z are coordinates of Ver1. x,y,z, respectively, stands for the particle Xi on x,y,z axis coordinate. x_step and z_step, respectively, stands for space interval on x,z axis coordinate of the terrain elevation. m_height stands for the two-dimension array which stores the terrain data:
(4.1)m=int_min(xix_step),n=int_min(ziz_step),y=m_height[m][n],Ver1(x,y,z)={m*x_step,y,n*z_step}.Step 2.
Judge whether the particleXi is in the area that rebounded by the terrain. First, calculate the spatial distance |Xi-Ti| from particles Xi to Ti. If |Xi-Ti|>R, then the particle Xi is not in rebound area, so go to Step 1, and deal with particle Xi+1. Otherwise, |Xi-Ti|≤R, then enter Step 3, R values for h which is the radius of the smoothed particle.Step 3.
Calculate force of rebound boundary. When fluid particleXi is close to terrain, the particle Xi will be rebounded by the terrain surface. The rebound force fi can be calculated by (4.2), where ui is speed of the particle, Kstiff is the hardness parameters of the terrain, Kdamp is deceleration parameter of terrain to the particle, and n is the normal unit vector of the triangular patches that include the intersection Ti. When the fluid particles gradually approach the terrain surface, the rebound force increases gradually until the velocity of the particle decreases to zero. At this time, the particle is still subjected to the effect of rebound which increases the velocity of the particle, thus preventing fluid particles to penetrate through the terrain. The direction of the rebound force is normal direction of the triangular mesh plane. Note that the rendering system only renders the clockwise triangular patches; we adopt clockwise triangles for the calculation of the normal unit vector of triangular patches that the particles locate in:
(4.2)fi={Kstiff*(2R-d)-ui*n*Kdamp*n,R>d,d=|Xi-Ti|,0,otherwise.


## 5. Visualization of the Fluid Surface

The fluid surface reconstruction is important for fluid simulation. Lighting and texture rendering of the fluid surface will improve the fidelity of fluid scene. Iwasaki et al. [15] proposed Point Splatting method to build the surface rendering, and this method has highly real-time rendering efficiency, but emptiness appears where fluid particles are relatively sparse [16–20]. In this paper, Marching Cubes algorithm is proposed to be applied in constructing fluid surface of free moving, which has the advantages of simple operation and fast drawing. The process is as follows.Step 1.
The 3D space is evenly divided into spatial grids. Traverse all the particles in the fluid region, calculate minimum coordinates (Vermin) and maximum coordinates (Vermax) of the particles which are bounding box of fluid region, and divide the box into homogeneous spatial grids.Step 2.
Calculate the density of the vertices of each spatial grids. Then use the divided grids as the sampling points of fluid scenes to get fluid density field data. Calculate the density value of each grid of the 8 vertices by using (2.1) and SPH interpolation on particles of each grid.Step 3.
Determine the threshold value of density isosurface of the fluid surface. According to the fluid pressure equal to atmospheric pressure, and the state equation of an ideal atmosphere, threshold value of density isosurface of the fluid surface can be obtained as follows:(5.1)ρsurface=pk,
where ρsurface stands for the fluid surface density, p for atmospheric pressure, and k for gaseous constant.Step 4.
Compare each vertex density value of spatial grid with threshold value of isosurface. If the vertex density value is greater or equal to the isosurface value, vertex value is 1 and the vertex is in the isosurface; otherwise is 0 and the vertex is outside of isosurface. According to the result of the comparison, structure the grid state table.Step 5.
Calculate the vertex of density isosurface coordinates of each spatial grid. If a vertex of a side of a grid is in the isosurface, while the other vertex is out of the isosurface, then, this side inevitably intersects with the desired isosurface. Firstly, according to the grid state table, get the grid’s sides which intersect with the isosurface. Then calculate the intersection of grid’s side and isosurface by linear interpolation method.Step 6.
Draw density isosurface. By using the central difference method and the linear interpolation method, calculate vertex normal of each triangle. Finally according to the coordinate values of each triangle vertex and normal vector, draw density isosurface.


## 6. Experiment Results

In this paper, the experiment platform is Intel Core 2 Duo CPU T6670@2.2GHZ. Main memory is 2 G, and graphic card is NVIDIA GeForce 9300GS with memory of 256 M, and operating system is Windows XP, and the drawing SDK is OSG (OpenSceneGraph).Figure3 shows fluid outflowing from a topographic position in 3DGIS, and its dynamic process of evolution on the terrain surface. The fluid particle number is 30000. Figure 4 shows the Figure 3 Scene by Marching Cube algorithm to construct the effect of fluid surface of free moving.Figure 3
Fluid in the dynamic evolution process of particle simulation (particle number 30000).Figure 4
Effect of reconstruction of fluid surface mesh reconstruction (particle number 30000).As this algorithm applies Marching Cube algorithm which is the same as [8] to construct the fluid surface and field particle search and complex boundary interaction algorithm which are used in this algorithm are different to [12], thus they are comparable. Table 1 shows the comparison of simulation speed of the algorithm and two other kinds of algorithm. Figure 5 shows comparison of computing time per frame. Table 1 and Figure 5 show that this algorithm reduces the computation time per frame.Table 1
Comparison of the three algorithms on simulation rate.
This paper algorithm
[12] Algorithm
[8] Algorithm
Particle number
FPS
Particle number
FPS
Particle number
FPS
2000
153.4
3000
68
4000
33
4000
94.3
8000
45.2
8000
26
16000
26.4
20000
12.7
20000
10
32000
14.5
64000
8.2
30000
6.7
128000
3.8Figure 5
Comparison of the three algorithms on simulation rate.


## 7. Conclusion

In order to simulate the fluid scene in 3DGIS, this paper proposed an SPH algorithm for fluid scene simulation. The experiment results show that this algorithm can calculate in real-time and has a good real-time rate of rendering frame, and achieve the desired goals. In near future we intend to further improve the algorithms of extracting fluid surface and enhance reconstruction speed and accuracy of the fluid surface mesh.



---

*Source: 635631-2012-10-22.xml*