Improvement in the Boundary Conditions of Smoothed Particle Hydrodynamics
T. Harada, S. Koshizuka and Y. Kawaguchi
The
University of Tokyo, Japan
takahiroharada@iii.u-tokyo.ac.jp
Contents
1. Introduction2. Related Work
3. Smoothed Particle Hydrodynamics
3.1. Governing Equations
3.2. Basics of Smoothed Particle Hydrodynamics
3.3. Discretization
3.4. Kernels
3.5. Neighboring Particle Search
4. Wall Weight Functions
4.1. Density Computation
4.2. Viscosity Term
4.3. Pressure Term
4.4. Computation of Wall Weight Functions
5. Results
6. Conclusions
References
Absrtact
In this paper, we propose an improvement in boundary condition computation of Smoothed Particle Hydrodynamics. Generally, to calculate a wall boundary in particle methods, a wall is converted into particle representation and used as static particles. However, it has several issues, e.g., a boundary is discretized into particles and thus does not result in accurate representation of the boundary shape. For precise representation of the shape, the resolution of a simulation must be increased and thus a large number of particles are required. The proposed method does not use boundary particles but rather uses a distance function to represent the boundary. Therefore, the shape of a boundary can be represented precisely without increasing the resolution of the particles. Thus, the total particle number can be reduced and it lighten computational burden. Hence, the proposed method is useful for real-time applications.
We come across with many kinds of fluids that have a complex behavior. Therefore it is difficult to manually create an animation of liquid motion. However, we can make such an animation by solving the governing equations. Therefore, fluid simulation is an important research theme in computer graphics. There are two major approaches to simulate fluid motion: Eulerian and Lagrangian methods. Eulerian methods use a grid to store physical values and discretize the governing equations on the grid. In contrast, Lagrangian method does not use a grid but uses particles, which does not have any fixed connectivity information between them and moves freely in the computation region. Physical values are stored on the particles.
Particle methods are well suited for a simulation of a free surface because particles themselves represent a fluid and they do not have to track the surface. Another benefit is that they are not accompanied by numerical dissipation caused by the convection term because the convection terms are calculated by moving the particles. When Eulerian methods are used, we have to overcome these issues.
Moving Particle Semi-implicit (MPS) method and Smoothed Particle Hydrodynamics (SPH) are particle methods used for fluid simulation[1][2]. The MPS method that can solve incompressible flow is well studied in the engineering field[3][4]. SPH, that computes compressible flow, is often used in computer graphics because computation cost of SPH is lower than that of the MPS method[5][6]. By lowering the compressibility, SPH can solve incompressible flow approximately.
In particle methods, wall boundaries are converted to particles and are used as boundary[7]. To prevent penetration of fluid particles through the wall boundary, several layers of boundary particles have to be prepared. The number of boundary particles is proportional to the surface area of the boundary. Therefore, the number of boundary particles, as well as the computational cost, increases with the increase in the complexity of the boundary shape. However, there is a drawback of using boundary particles. The shape of the boundary cannot be represented precisely. Therefore it is difficult to compute fluid motion in complex shapes with boundary particles.
In this paper, we present an improved computation model for boundaries, which does not use boundary particles but uses a distance function calculated from a polygon model. Wall weight functions are introduced to compute physical values of fluid particles close to a boundary. The benefit of the present method is that the shape of a boundary is precisely represented which is difficult when boundary particles are used. As a result, a fluid motion in complex shapes can be simulated easily. The proposed method reduces the computational burden because computation of physical values of boundary particles is no required and the wall weight functions can be precomputed.
This paper is
organized as follows. After the discussion about related work, the basics of
SPH and its formulations are described. Then, we present the wall weight
functions. Finally, several examples are presented in Section 5 and conclude
the paper with Section 6.
SPH was developed
and have been studied for astrophysical applications to simulate astral bodies[8]. Monagham et al. applied SPH to fluid
simulations[2].
Several kernels and formulations of SPH have been developed in engineering
fields[9].
Morris et al. studied multi-phase fluid simulation with surface tension[10]. However, SPH cannot solve incompressible
fluid because it does not solve the mass conservation equation. SPH solves a
fluid by lowering the compressibility. On the other hand, the MPS method
composes a Poisson equation with particle density and solves it for the
pressure[1].
Therefore, the MPS method can be used for incompressible fluid computation,
although it is much expensive than SPH. Later, several researchers employed the
similar approach to the MPS so that SPH could be used for incompressible flow[11][12].
In computer
graphics field, Desbrun et al. used SPH to simulate the motion of soft
substances[13].
Muller et al. used SPH for a real-time simulation of a liquid because it has
several advantages over the Eulerian methods such that it is easy to compute
free surfaces[5].
They used several-thousands of particles and then coupled with mesh-based
elastic body simulation and multiphase flow[14][15]. Premoze et al. introduced the MPS method
to the graphics community[6].
To improve the
performance of a simulation, several researchers used Graphics Processing Units
(GPUs) to accelerate computations. Kolb et al. accelerated particle simulation
without the interaction between particles[16]. They showed that GPUs have the potential to drastically accelerate a
particle-based simulation. Kipfer et al. sorted particles by their coordinates
and solved collision between them[17]. However, their method cannot deal with perfect collisions. Harada et al.
developed a technique to perform all the computations of SPH using a GPU, and
they simulated about 60K particles in real-time[18]. In addition, they also studied Sliced
Data Structure which improved the memory efficiency of the grid used for the
acceleration of neighboring particle search[19].
3.
Smoothed Particle Hydrodynamics
A
fluid is incompressible, so the pressure force works as a force to make the
fluid incompressible. The viscosity force makes the velocity of a fluid
uniform. The variation of velocity of a fluid is equal to the force on the
fluid. Therefore, the momentum conservation equation is
(1)
where 
are the density, velocity, pressure, viscosity
coefficient and external force, respectively. The terms on the right side of
the Eqn. (1) are pressure, viscosity and external
force terms, respectively. The term on the left side is the Lagrangian
derivative, which can be expressed as
. (2)
The second term in Eqn. (2) is the advection term. When a grid is used to compute the advection term, the physical values at any point have to be interpolated from the grid values, which results in a numerical diffusion. However, the Lagrangian derivative is a derivative on a flow, i.e., on a moving particle: therefore, particle methods can evaluate it by advecting particles. This means there is no numerical dissipation in the advection computation.
3.2. Basics of Smoothed Particle Hydrodynamics
SPH is a numerical method that discretizes a fluid into particles. To solve Eqn. (1), a physical value
at an arbitrary position must be calculated. In SPH, a physical value 
at position 
is obtained by integrating
physical values over the space as follows:
(3)
where W is a normalized kernel. Therefore, the integration of the kernel W over the space is unity.
In SPH, particles store physical values and
it is distributed around each particle as shown in Figure 1. Therefore, the
integral in Eqn. (3), which is the value at , is the summation of the values 
on the neighboring particles. Therefore, Eqn. (3) can be expressed as
(4)
where 
are the mass and density of particle j,
respectively and 
. In general, kernel W is
defined as a function that is equal to zero at positions far from the particle. Therefore, to calculate a physical value 
of a particle i,
physical values at particles close to particle i are required.
For example, the density at can be calculated as follows:
. (5)
We are now able to calculate a physical
value at an arbitrary position. However, the gradients of physical values on
particles also have to be calculated. The gradient 
of at is calculated by taking the
gradient of Eqn. (3).
(6)
Here, the partial integral is used and the term to the right side can be decomposed to two terms as follows.
(7)
The first term on the right side can be zero using Gauss's low. Therefore, the gradient of a physical value is calculated as
. (8)
When we use particles, the expression can be converted as follows.
(9)
Therefore,
to calculate the gradient of a physical value, we have to differentiate the
kernel and not the physical values.
Using Eqn. (9), the pressure term can be expressed as
(10)
where 
is the pressure on particle j. However, this equation makes the force between
particles i and j asymmetrical. Therefore, the pressure force is modeled with the pressure of
particles i and j.
(11)
Although there are several approaches to
compute the pressure, we adopted a simple one. From the state equation of the
gas, the pressure at particle i can be
calculated as
(12)
where 
is the rest density which is the density at
the initial particle distribution.
To compute the viscosity term, laplacian operator has to be defined for the particles. Here, we first describe the formulation of the laplacian operator.
A physical value 
at particle i
is approximated with the physical value at particle j by
(13)
Therefore, 
at particle j is
. (14)
The laplacian is
the divergence of the gradient (
).The divergence operator can be
modeled as Eqn. (9) and substituted it into Eqn. (14).
(15)
Here, we define the kernel of the viscosity 
, and Eqn
(15) is expressed as follows.
(16)
Using Eqn. (16), the viscosity term can be calculated as
. (17)
Kernels are required to calculate several
terms. First, the kernel used for density computation is defined as
(18)
where d is the effective radius. Neighboring
particles in the effective radius of a particle have an effect on it but the
particles outside of the effective radius have no effect on the particle.
Although the basics of SPH state that the differentiated values are calculated
by differentiating the kernel, the gradient of Eqn.
(11) decreases as the particle approaches. This means that the pressure force
decreases as the particle comes closer. In reality, the pressure has to
increases as the particle approaches. Therefore, we used the following kernel
for the pressure term as used by Muller et al.[5].
. (19)
For the viscosity term, we used this kernel as follows.
(20)
3.5. Neighboring Particle Search
To compute a physical value stored on a
particle, the neighboring particles have to be searched for. In particle-based
methods, particles move freely and therefore the neighboring particles have to
be searched for in every time step. If the brute-force strategy is employed,
the computation would be 
, which is not practical for a simulation with a large number of
particles. Therefore, we introduced a uniform grid to make the search
efficient. Before the search, a particle index is stored in the voxel to which the particle belongs. With this grid, the
search becomes 
because the indices of the neighboring
particles are stored in the neighboring voxels.
In this study, we used a uniform grid:
therefore, the sliced data structure proposed in [19] can be
used.
When simulating a fluid, boundary conditions have to be considered. One choice is to generate boundary particles and treat them as static particles. If we employ this strategy, several layers of boundary particles have to be prepared to precisely compute the density and pressure of the boundary particle. Therefore, a large number of particles are required other than fluid particles. Furthermore, boundary particles cannot precisely represent the shape of a boundary. Even if the boundary is flat, it will be considered as a bumpy boundary. When we simulate a fluid on a slope, the fluid does not flow smoothly.
The improved model for boundary conditions
overcomes these issues. The method introduces wall weight functions and
computes the contribution of a boundary using the distance to the boundary.
Therefore, a fluid flows smoothly on a slope because the effect of the boundary
to points, which have the same distance to the boundary, is computed as the
same. What affects the fluid particles are the density, viscosity term and
pressure terms. We explain these models in the following subsections.
The density of a fluid is calculated as Eqn. (5) when we use wall particles as the wall boundary. We assume that there is a wall boundary within the effective radius of particle i. The computation of fluid density is divided into the contribution of the fluid and wall particles as follows.
(21)
The first and second terms on the right
side are the contribution of fluid particles and wall particles, respectively.
To compute the second term, temporary wall particles are placed perpendicular
to the perpendicular line from particle i to the wall boundary as shown
in Figure 2. The distribution of wall particles in the effective radius is
determined uniquely when the vector from the wall 
is determined.
Since the distance 
is used in the
calculation of the weight function W, the contribution of wall particles is calculated as a function of
the distance to the wall boundary as follows.
(22)
We refer to 
as the wall weight function for density.
The viscosity term is also divided into the
contribution of fluid particles and that of wall particles. The viscosity force 
for particle i can be divided into the
viscosity force from fluid particles 
and that from wall particle 
as follows.
(23)
We
assume that the velocities of all the wall particles are the same. The
viscosity force from the wall is calculated as
(24)
We
are dealing with an incompressible fluid. Although the mass conservation
equation is not solved in SPH, it can be used for near-incompressible fluid
computation. This means that the density inside of the fluid is approximately
uniform. When wall particles are used for the boundary conditions, the density
of wall particles is also approximately uniform. Therefore, in Eqn. (24),
of all the wall
particles is constant by assuming that the density of wall particle 
is constant. The
distribution of wall particles in the effective radius is also determined
uniquely when
is determined.
Therefore, the contribution of wall particles in the viscosity term is
calculated as follows by using the distance to the boundary 
from particle i.
(25)
We refer to 
as the wall weight function for viscosity.
When an incompressible fluid is solved by
using a particle method, the pressure term works to make the particle number
density constant. Therefore, the pressure force from wall 
to fluid particle i is modeled as the force that
makes the distance between the fluid particle and the wall boundary as the rest
distance d.
When particle i is located at a the position
where the distance to the wall boundary is (is smaller than d), the pressure force pushes the particle i at a distance 
in the direction of 
, which is the normal vector of the closest wall
boundary point at 
.
Thus the pressure force from the wall is modeled as
. (26)
Fig. 2: Wall weight functions
4.4. Computation of Wall Weight Functions
Since the wall weight functions in Eqns. (22) and (25) depends on the distance to the wall boundary , all we have to calculate the distance in order to compute
these wall weight functions. The wall weight functions can be computed in advance
and the values are referred at the computation of the fluid. Particles are
placed as shown in Figure ¥ref{img:weightFunction}
and the wall weight functions are calculated by integrating the weight
functions numerically because they cannot be computed analytically. The wall
weight function is computed at a few points within the effective radius. At the
other points, the function is interpolated from the calculated values.
By precomputing
the wall weight functions, the contribution of the wall is computed easily by
referring the value corresponding to the distance. However, to calculate the
distance to the boundary, we have to compute the distance to all the polygons
comprising the boundary and select the minimum value. When the number of
polygons increases, the computation becomes too expensive for real-time
applications. Therefore a distance function is introduced to reduce the
computational burden. As the distance function is the distance to the closest
boundary point, the distance between particle i and the
boundary is obtained from the
distance function by using particle position 
. When the wall boundary does not deform, the distance
function does not change. Therefore, the distance function can also be computed
in advance. Consequently, we do not have to compute the distance to all the
polygons but we refer the precomputed distance
function at the computation of fluid. When the pressure term is calculated, the
normal vector of the boundary point which is the closest to particle i is required as well as the distance to the boundary. The vector
from the closest point on the boundary has the same direction as the normal
vector,
which is called the gradient vector. If the gradient
vector is obtained, we can calculate the normal vector. In addition, the
gradient vector does not change if the wall boundary does not deform. Hence the
gradient vector can be calculated beforehand and utilized for fluid computation
as well. Also, the distance function is the absolute value of the gradient
vector because it is the vector to the closest point on the boundary.
Fig. 3: A computation result of free surface flow. A gargoyle model was used as a boundary
We implemented the proposed method on a PC with Pentium 4 3.6GHz CPU and 3.0GB RAM. The programs were written in C++. Several types of polygon models were used as wall boundaries and some results are shown in this paper. The gradient vector was precomputed and interpolated at run time. Since the wall weight functions were also precomputed, the computation cost of the proposed wall boundary calculation method is low. As we described above, the method can represent the shape of a boundary precisely without increasing the number of particles. To precisely represent the shape with wall particles, we have to increase the resolution of particles, which results in a large number of boundary particles.
Figure 3 is a result of a simulation of pouring fluid into a gargoyle model (upside down). Surface polygons were extracted from the particles used in the simulation using Marching Cubes[21] and were rendered. The model has two thin wings, therefore it is difficult to simulate in such regions with boundary particles. However, we can see that the fluid flowed into these regions without any artifacts. This is a benefit of the proposed method. In Figures 4 and 5, a liquid was poured into a Buddha and dragon models. Although the shapes of these two models are very complex, the fluid flows smoothly in these models. These results are rendered with matte shader in OpenGL in real-time. The computation times with these polygon models as wall boundaries are shown in Table 1. We can see that one time step is finished within a second.
When the wall particle boundary is used, the number of wall particles which have to be generated increases as the geometry of wall boundary becomes more complicated and the surface area increases. Wall particles are generated from the gargoyle, dragon and Buddha models by using surface voxelization method[20]. Table \2 shows the number of wall particles, total number of particles and their ratio. We can see that a large number of particles were required when boundary particles were used. Since the proposed method does not require boundary particles, nearly half of the particles are eliminated.
We also implemented an SPH simulation entirely on the GPU based on the study of Harada et al.[18] and implemented the boundary
computation method. The proposed method is also useful for the simulation on
the GPU because the video memory on a GPU is smaller than the main memory and thus we can use more fluid particles than
a simulation with boundary particles. In Figure 6, droplets are thrown into a tank. Total number of
particles is about 1M. Figure 7 shows
the result of pouring a fluid onto another type of fluid.
Figure 6 is
rendered with pointsprites in OpenGL and Figure 7 is
rendered after the simulation with a ray tracer.
The proposed method still has certain limitations. The method assumes that all
the boundaries are flat and does not consider the curvature of the boundary.
Although we can obtain plausible results and observe no artifact from this, the curvature should be considered for much accurate computations.
|
Model |
Fluid Particles |
Computation
Time |
|
Gargoyle |
20,000 |
309.4 |
|
Dragon |
20,000 |
281.2 |
|
Buddha |
20,000 |
296.8 |
|
Model |
Wall |
Total |
Ratio |
|
Gargoyle |
26,688 |
46,688 |
0.572 |
|
Dragon |
18,582 |
38,582 |
0.482 |
|
Buddha |
12,084 |
32,084 |
0.377 |
Fig. 4: A computational result of free surface flow in a Buddha model
Fig. 5: A computational result of a free surface flow in a dragon model
In this study, we presented an improvement in boundary condition computation in SPH. The method uses distance function instead of boundary particles. This simplifies the simulation and enables us to simulate a fluid motion in complex shapes with ease. We obtained several results and proved that the proposed method can be used to reduce the total number of particles.
In this study, the curvature of the
boundary is not considered
and is a topic for future studies. Another possibility
would be to
extend our study for a moving boundary. Interaction between elastic and rigid boundaries would be challenging.
Fig. 7: A result of two phase flow. The tank is filled with a fluid and another fluid is poured in the tank
[1] S. Koshizuka and Y. Oka,
Moving-Particle Semi-Implicit Method for Fragmentationf
of Incompressible Flow, Nucl. Sci.
[2] J.J. Monaghan, Simulating Free Surface Flows with SPH, Journal of Computational Physics, Vol. 110, No. 3, pp. 333-353, 2000
[3] K. Shibata,
[4] T. Harada, Y. Suzuki, S. Koshizuka, T. Arakawa and S. Shoji, Simulation of Droplet Generation in Micro Flow using MPS Method, JSME International Journal Series B, Vol. 49, No. 3, pp. 731-736, 2006
[5] M. Muller, D. Charypar, and M. Gross, Particle Based Fluid Simulation for Interactive Applications, Proc. of Siggraph Symposium on Computer Animation, pp. 154-159, 2003
[6] S. Premoze, T. Tasdizen, J. Bigler, A. Lefohn and R. Whitaker, Particle-based Simulation of Fluids, Computer Graphics Forum, Vol. 22, No. 3, pp. 401-410, 2003
[7] J. Morris, P. Fox and Y. Zhu, Modeling Low Reynolds Number Incompressible Flows using SPH. Journal of Computational Physics, Vol. 136, pp. 214-226, 1997
[8] L.B. Lucy, Numerical Approach to Testing The Fission Hypothesis, Astronomical Journal, Vol. 82, pp.1013-1024, 1997
[9] J.J. Monaghan, Smoothed Particle Hydrodynamics, Annu. Rev. Astrophys., Vol. 30, pp.543-574, 1992
[10] J.P. Morris, Simulating Surface Tension with Smoothed Particle Hydrodynamics, International Journal of Numerical Methods in Fluids, Vol. 33, No. 3, pp. 333-353, 2000
[11] S.J. Cummins and M. Rudman, An SPH Projection Method, Journal of Computational Physics, Vol. 152, No. 2, pp. 584-607, 1999
[12] Y.M.L. Edmond and S. Songdong, Simulation of Near-shore Solitary Wave Mechanics by an Incompressible SPH Method, Applied Ocean Research, Vol. 24, pp. 275-286, 2002
[13] M. Desbrun and M. Gascuel. Smoothed Particles: A New Paradigm for Animating Highly Deformable Bodies, Proc. of Eurographics Workshop on Computer Animation And Simulation, pp. 61-76, 1996
[14] M. Muller,
[15] M. Muller, B. Solenthaler, R. Keiser and M. Gross, Particle-based Fluid-fluid Interaction, Proc. of Siggraph Symposium on Computer Animation, pp. 237-244, 2005
[16] A. Kolb, L. Latta and C. Rezk-Salama, Hardware Based Simulation and Collision Detection for Large Particle Systems, Proc. of the ACM Siggraph/ Eurographics Conference on Graphics Hardware, pp. 123-131, 2004
[17] P. Kipfer, M. Segal and R. Westermann, Uberflow: A GPU-based Particle Engine, Proc. of the ACM Siggraph/Eurographics Symposium on Computer Animation, pp. 23-31, 2002
[18] T. Harada, S. Koshizuka and Y. Kawaguchi, Smoothed Particle Hydrodynamics on GPUs, Proc. of Computer Graphics International, pp. 63-70, 2007
[19] T. Harada,
[20] Z. Dong, W. Chen, H. Bao, H. Zhang and Q. Peng, Real-time Voxelization For Complex Polygonal Models, Proc. of Pacific Graphics, pp. 73-78, 2004
[21] W. Lorensen, H. Cline, Marching Cubes: A High Resolution 3D Surface Construction Algorithm, Proc. of the 14th Annual Conference on Computer Graphics and Interactive Techniques, pp. 163-169