Fast Simulation of Viscous Fluids with Elasticity and Thermal Conductivity
Using Position-Based Dynamics
Tetsuya Takahashia,b , Tomoyuki Nishitab,c , Issei Fujishiroaa Keio University
b UEI Researchc Hiroshima Shudo University

Abstract
Viscous fluids are ubiquitous, and reproducing their damped motions has been in demand for many applications. The most prevalentapproach to simulating viscous fluids is based on the Navier-Stokes equations and necessitates viscosity integration. However, tosimulate viscous fluids in a numerically stable manner, using explicit viscosity integration severely restricts time steps and requires
an excessively long period for computation. In this paper, we propose a novel particle-based Lagrangian method for eﬃciently
simulating viscous fluids by adopting position-based constraints. Our method uses the geometric configuration of particles for the
positional constraints to approximate the dynamics of viscous fluids using position-based dynamics; thus the method can plausibly
generate their motions while allowing for the use of much larger time steps than those previously adopted in the viscous fluid
simulations. We also propose an associated boundary-handling scheme for position-based fluids to precisely specify boundary
conditions for the constraints. Additionally, we reproduce elastic deformations of materials by controlling the constraints andincorporate thermal conduction into our framework to simulate resultant changes in particle properties and phase transition in thematerials. By adjusting parameters, our method can encompass complex motions of fluids with diﬀerent properties in a unified
framework. Several examples demonstrate the eﬀectiveness as well as versatility of our method.
Keywords: Fluid simulation, viscous fluid, position-based dynamics, geometric constraint, elasticity, thermal conductivity.

1

1. Introduction

dominant factor for determining time steps in the viscous fluid
simulations; as a result, ∆t ≤ O(ρ∆x2 /µ) is necessary at the
29 significant sacrifice of computational eﬃciency.
30
In this paper, therefore, we propose a novel particle-based
31 Lagrangian method for simulating viscous fluids faster by al32 lowing for the use of larger time steps with position-based con33 straints. The core idea of our method is to take full advantage of
34 the geometric configuration of particles for the positional con35 straints to approximate the dynamics of viscous fluids. Since
36 our method does not require solving the diﬀusion equations,2
37 we can ignore the restriction ∆t ≤ O(ρ∆x /µ); thus, time steps
38 with our method can be much larger than the limited time steps
39 for explicit viscosity integration.
40
We use position-based fluids of Macklin and Müller [8] as
41 our fluid solver, where fluid is simulated using position-based
42 dynamics [9, 10], which directly corrects particle positions to
43 produce plausible motions of objects. To precisely specify bound44 ary conditions with the positional constraints, we also propose
45 an associated boundary-handling scheme for position-based flu46 ids with non-fluid particles on object surfaces. Although the
47 position-based method is less accurate than the traditional force48 based one, our method can generate visually plausible viscous
49 fluid motions. Figure 1 illustrates the eﬀect of our method with
50 an example of melted chocolate.
51
As an extension to [11], we newly added handling of elas52 tic materials and thermal conduction to our framework and im53 proved the constraint for viscosity to make it easier to produce27

28

Viscous fluids are often found in our daily lives, and theircomplicated behavior has fascinated us. Unlike inviscid fluids,
4 viscous ones, such as honey, melted chocolate, wax, paints, sap,
5 slime, lava, and blood, exhibit characteristic and interesting mo6 tions. Since their damped motions produce visually attractive
7 phenomena and play important roles in enhancing visual real8 ity, these viscous fluids have frequently been adopted in movies
9 and video games.
10
The most prevalent approach to simulating viscous fluids
11 is based on the Navier-Stokes equations and has successfully
12 simulated viscous fluids [1, 2, 3, 4, 5].
In computer graph13 ics, the equations are usually solved with three diﬀerent ap14 proaches: Eulerian, Lagrangian, or Eulerian-Lagrangian hy15 brid.
We herein focus on the particle-based Lagrangian ap16 proach due to its versatility.
17
To solve the Navier-Stokes equations with particle-based
18 methods, such as Smoothed Particle Hydrodynamics (SPH) [6],19 explicit viscosity integration is often used due to its simplic20 ity and eﬀectiveness. However, this explicit scheme is likely to
21 fail in plausibly reproducing highly viscous fluids under high22 spatial resolution because, to describe viscous eﬀects based on
23 diﬀusion equations, we need to satisfy the numerical stability
2
24 condition µ∆t/(ρ∆x ) ≤ 1/2 (µ: dynamic viscosity, ∆t: time
25 step, ρ: fluid density, and ∆x: interparticle distance), accord26 ing to von Neumann stability analysis [7]. This condition is a2
3

Preprint submitted to Computers & Graphics

July 7, 2014

Figure 1: Melted chocolate running down the surface of a white ball. With our method, chocolate flows on and sticks to the ball, and slowly flows over the floor,
due to the geometric constraints between particles, which approximate the dynamics of viscous fluids by restricting movements of the particles.

54
55

56
57
58
59

60
61
62

63
64
65

66

desirable fluid motions by reducing the number of parameters.
In summary, our key contributions are three-fold:• Much larger time steps are available to simulate viscous
fluids than those previously used to numerically stabilize the simulation with explicit viscosity integration by
avoiding the strict time step restriction ∆t ≤ O (ρ∆x2 /µ);• The associated boundary-handling scheme for positionbased fluids properly addresses the position-based constraints on object surfaces; and• Various behaviors of fluids with diﬀerent properties dueto thermal conduction can be generated and controlled byadjusting parameters in a unified framework.
2. Related Work

We briefly explain methods for simulating viscous fluids
and techniques for particle-based methods.
69
Viscous fluids As a versatile framework for various types of
70 materials, the Lagrangian finite element method has been pro71 posed and accurately simulated the dynamics of viscous fluids
72 [12, 13, 14, 15]. However, this method is inappropriate for fast
73 simulation due to the high computational cost of remeshing and74 solving linear systems.
75
Bergou et al. [16] and Batty et al. [17] proposed a dis76 crete model specialized for describing dimensionally reduced
77 materials and successfully simulated the distinctive behavior of
78 viscous threads and sheets, respectively, yet their methods are
79 inapplicable to three-dimensional viscous volumes.
80
Desbrun and Gascuel [18] adopted SPH and simulated vis81 cous materials with an artificial viscosity term [19]. This term
82 was also used by Stora et al. [20] for lava simulation. Müller
83 et al. [6] proposed a viscosity term denoted by the Laplacian
84 operator and simulated dynamic fluid with viscosity.
Solen85 thaler et al. [21] combined this term with an elastic force term
86 to simulate fluid with viscosity and elasticity and simulated
87 changing particle properties and phase transitions in materi88 als by solving heat equations in a unified framework. Becker
89 et al. [22] improved Solenthaler et al.’s method [21] to han90 dle rotational motions of elastic objects. Paiva et al. [5] used a
91 generalized Newtonian model for viscous fluid simulation and
92 adopted XSPH [23], which smoothes velocities for coherent
93 particle movements. For simulating viscoelastic materials, Mao
94 and Yang [24] introduced a new elastic force term, which is a67
68

nonlinear corotational Maxwell model, into the Navier-Stokes
equations. Chang et al. [25] used a simplified Maxwell model,
97 which drops a rotational tensor, and simulated viscoelastic ma98 terials including the eﬀect of thermal conduction. A Maxwell
99 model was also used to simulate the coiling phenomenon [26].
100
Recently, methods that combine SPH with another tech101 nique have been proposed.
Gerszewski et al. [27] presented
102 a method that uses an aﬃne transform to approximate the mo103 tions of neighboring particles for reproducing elastoplastic ma104 terials. Takamatsu and Kanai [28] simulated viscoelastic ma105 terials by combining SPH with shape matching. Dagenais et
106 al. [29] reproduced viscous fluid motions in SPH by adding ex107 tra forces that move particles to their original positions.
108
Several intuitive methods have also been proposed for vis109 cous materials. Miller and Pearce [30] and Terzopoulos et al. [31]110 reproduced viscous materials using a spring-based model that111 computes the repulsion and attraction forces between particles.
112 This idea was also adopted by Steele et al. [32] and Tamura
113 et al. [33]. Clavet et al. [34] extended this model to simulate
114 materials that exhibit elasticity, viscosity, and plasticity.
115
Since particles are controlled on the basis of their geomet116 ric relations without estimating force fields, our method can be
117 categorized as a spring-based method. However, our method
118 diﬀers from the spring-based method in that ours is position119 based; thus ours can perform numerically stable simulation with120 larger time steps.
121
Particle-based methods The particle-based method for sim122 ulating fluid is popular in computer graphics applications and123 has been developed for enforcing fluid incompressibility, one124 way and two-way solid coupling, and fluid-fluid interactions
125 [35]. For incompressible fluids, various pressure solvers have
126 been proposed, e.g., Weakly Compressible SPH [36], Predictive127 Corrective Incompressible SPH [37], and Implicit Incompress128 ible SPH [38]. Interactions between fluid and objects are im129 portant and have also been improved for one-way and two-way
130 coupling with rigid bodies [39, 40, 41] and deformable objects
131 [42].
Additionally, for fluid-fluid interactions, several meth132 ods have been proposed to generate bubbles [43], simulate fluid
133 with significant diﬀerences in density [44], and produce the ef134 fects of diﬀusing, cleansing, and foaming [45].
135
To accelerate particle-based fluid simulations, sophisticated
136 techniques have been presented for multi-core CPUs [46] or
137 GPUs [47, 48]. As a diﬀerent direction to eﬃcient computation,
138 adaptively sampled particles reduced computational cost while
139 preserving the quality of the simulation results [49, 50, 51].
95
96

2

where ρ0 is the rest density of the fluid, and ρi is estimated withthe summation approach commonly used in SPH, taking non141
We use position-based fluids [8] as our underlying fluid 184 fluid particles into account [41]:142 solver because we can adopt larger time steps than those used∑
∑ ρ0
143 in SPH; thus, we can take full advantage of our geometric apρi =m j Wi j +
Wik ,
(5)
δk
144 proach to maximize time steps.
jk140

3. Proposed Method

182
183

We briefly summarize the fundamental formulations of position185 where j denotes a neighboring fluid particle of i, k a neighborbased dynamics [9, 10] in Section 3.1 and explain our fluid
186 ing non-fluid particle, Wi j the short for a kernel W(pi − p j , h)∑
147 solver in Section 3.2. In Section 3.3, we give details on how to
187 with a kernel radius h, and the number density δi =k Wik . The
148 simulate viscosity in our framework. We explain the handling
188 gradient of constraint Eq. (4) with respect to particle l is given
149 of elasticity in Section 3.4 and give an algorithm for controlling
189 by150 constraints in Section 3.5. We explain thermal conduction and1 ∑∑ 1
151 its eﬀects in Section 3.6 and show the procedure of our method
if l = i,j m j ∇Wi j + k δk ∇Wikρ0

152 in Section 3.7.
 mj
∇pl Cdens,i (p) = − ρ0 ∇Wi j
if l = j, (6)




1

−
∇W
if l = k,ik
153 3.1. Position-Based Dynamics
δk154
In position-based dynamics, an object is represented by a 190 and a density scaling factor λdens,i is computed by
155 set of N particles and a set of scalar constraints.
A particle
wiCdens,i (p)156 i has its mass mi , position p , and velocity vi .
The positioniλdens,i = − ∑
.
(7)
2
w
157 of each particle is iteratively corrected by resolving the conl l ||∇pl C dens,i (p)||158 straints, conserving the linear and angular momenta of the ob191 To handle the particle clustering induced by the negative pres159 ject. Let p denote the collection of particle positions p1 , . . . , p N ,192 sures around free surfaces due to the deficiency of particles,
160 and a constraint of p as C(p), which is preserved to satisfy the
193 Macklin and Müller [8] introduced an artificial pressure term:
161 bilateral (C(p) = 0) or unilateral (C(p) ≥ 0) condition.
()
162
We aim to find a correction vector ∆p such that C(p + ∆p) =Wi j 42
scorr,i = −k s h,(8)
163 0. This constraint is approximated byW(q, h)C(p + ∆p) ≈ C(p) + ∇pC(p)∆p = 0,
(1) 194 where k s is a parameter to control the eﬀect of artificial pres145
146

sure, and q is a vector whose norm is less than h. In summary,with Eqs. (2), (6), (7), and (8), the density correction vector is
197 computed by∆p = λ∇pC(p).
(2)
1 ∑m j (λdens,i + λdens, j + scorr,i )∇Wi j
∆p
=
dens,
i
166 For a particle i, we obtain a scaling factor λi with Eqs. (1) andρ0 j
167 (2):∑ 1+(2λdens,i + scorr,i )∇Wik .
wiCi (p)δkλi = − ∑
,
k
2
w
||∇
C
(p)||
pj i
j j
198 Additionally, XSPH viscosity [23] is used for coherent motions
168 where wi = 1/mi , and a position correction vector is given by
199 of particles with vi j = vi − v j :wiCi (p)∑ mj∑ 1∆pi = − ∑∇ C (p).
(3)
2 pi ivi ← vi + ϵav ji Wi j + ϵb
vki Wik ,j w j ||∇p j C i (p)||
ρj
δk
j
k
169
When a constraint with the unilateral condition is true, po170 sition correction is simply skipped. Note that the dimensionali- 200 where ϵa and ϵb are parameters for viscous eﬀects.
201
Since our aim is to simulate viscous fluids, we do not use
171 ties of C(p) and λ vary depending on the types of constraints.
202 vorticity confinement unlike in the original work [8]. Note that172 3.2. Position-Based Fluids with Boundary-Handling
203 by taking non-fluid boundary particles into the density estima173
Since the original method of position-based fluids [8] fails 204 tion in Eq. (5), they are also incorporated into the relevant for174 to correctly estimate particle density on object surfaces due to 205 mulations without any need of special treatments; therefore our
175 the deficiency of particles, we improve density estimation with 206 method ensures appropriate handling of particles on object sur176 non-fluid particles on surfaces by modifying the method of Ak- 207 faces.
208
Figure 2 compares the original position-based fluids adopt177 inci et al. [41] for position-based fluids.
178
In position-based fluids, constraints are imposed on each 209 ing level set boundary-handling [8] with our method. While
179 particle of a fluid. For a particle i with its density ρi , the dimen- 210 the original method [8] causes the artifacts of particle stack180 sionless density constraint with the bilateral condition is defined 211 ing along the vertical edge and unnatural particle aggregations
212 along the horizontal edges caused by erroneously underesti181 as
213 mated particle density, our method can generate natural fluid
ρi− 1 = 0,(4) 214 motions thanks to the proposed boundary-handling scheme.
Cdens,i (p) =ρ0and ∆p is restricted in the direction of ∇pC(p) with a scaling
165 factor λ to conserve the linear and angular momenta:

164

195
196

3

ǆƚĞŶƐŝŽŶ
௜௝

௜௝

(a) Position-based fluids [8](b) Our method

Figure 3: Illustration of extending interparticle distance di j for viscosity constraints. When two particles (green and orange) are connected, orange particle
can move within red area in left figure. After interparticle distance di j is updated with extension coeﬃcients βi and β j , orange particle can move within red
and blue areas in right figure.

Figure 2: Comparison of dam break simulation in a virtual box.

3.3. Viscosity
As described in Section 1, the current mainstream for sim217 ulating viscous fluids is based on the Navier-Stokes equations.
218 To describe viscosity in the equations, several viscosity terms
219 have been proposed [5, 6, 19]. However, their explicit schemes
220 fail to simulate viscous fluids under high spatial resolution in a
221 numerically stable manner due to the stability condition caused
222 by diﬀusion equations. To avoid this problem, we herein pro223 pose a particle-based method with geometric constraints.

௜  ௝
௜௝
2

215
216

3.3.1. Viscosity Constraint Mechanism
We approximate the eﬀect of viscosity using constraints be226 tween two particles connected to each other (hereafter, we call
227 these particles connected particles). Since particle velocities are
228 updated with position diﬀerences between consecutive simula229 tion steps in position-based dynamics, we can smooth particle
230 velocities by confining motions of connected particles within
231 certain distances, and this smoothing eﬀect is almost equivalent
232 to dissipating particle velocities to neighboring particles, as de233 scribed with the Laplacian operator. Specifically, we preserve
234 interparticle distances with constraints that directly correct par235 ticle positions for the smoothing eﬀect.
236
We dynamically generate, modify, and delete constraints
237 between particles in the simulation. If the distance of two par238 ticles is less than the kernel radius h, and there is no constraint
239 between them, we generate a new constraint. Since positions
240 of particles are corrected to keep their interparticle distance, we
241 modify a constraint by extending the distance so that particles
242 can move in the opposite direction from the other. We delete the
243 constraint between particles when the distance of the connected
244 particles is larger than H (we set H = 2h). After we explain
245 elasticity in Section 3.4, we give an algorithm for controlling
246 constraints of viscosity and elasticity in Section 3.5.
224

225

where pi j = pi − p j , and di j is initialized as the distance ofthe two particles when the constraint is generated. We adopt
255 the unilateral condition for the constraint because a constraint
256 with the bilateral condition tends to disturb the homogeneous
257 particle distributions induced by the artificial pressures. Vis258 cosity constraints handle only the expansion of fluid volumes,
259 and density constraints address the compression of the volumes.
260 Note that the dimensionality of the viscosity constraint diﬀers
261 from that of the density constraint. If Eq. (9) does not hold, we
262 correct the positions of connected particles using Eq. (3) with
263 a scaling factor sv (0 ≤ sv ≤ 1) to maintain their interparticle
264 distance:
pi j
sv wi
∆pvisc,i = −(||pi j || − di j ),(10)
wi + w j||pi j ||sv w jpi j∆pvisc, j = +(||pi j || − di j ).
(11)
wi + w j||pi j ||253
254

Since Eq. (9) is likely to be exceedingly strict for modifying
the viscosity constraint, we subtract the term αdi j to weaken
267 the condition as265
266

Dvisc (pi , p j ) = Cvisc (pi , p j ) − αdi j= di j − ||pi j || − αdi j ≥ 0.

(12)

If Eq. (12) does not hold, we extend an interparticle distanceto enable connected particles to slightly separate from the other270 (see Figure 3):268269

di j ← di j +βi + β j
di j ,
2(13)

where βi and β j are extension coeﬃcients that control the variation in the interparticle distance. Note that we modify the exten273 sion coeﬃcient, not viscosity value, and a high (low) extension
274 coeﬃcient is equivalent to low (high) viscosity. Our method can
247 3.3.2. Viscosity Constraint
275 produce motions of fluids with diﬀerent viscosity values using
248
To restrict particle motions when particles are separating
276 diﬀerent extension coeﬃcients, as shown in Figure 4.
249 from each other, we correct particle positions based on the vis250 cosity constraint C visc (pi , p j ) with the unilateral condition, de277 3.4. Elasticity
251 fined with the positions of two particles and their interparticle
278
Elasticity, which diﬀers from viscosity, is a characteristic
252 distance di j as
279 disposition that some materials exhibit. According to Zhou etCvisc (pi , p j ) = di j − ||pi j || ≥ 0,
(9)
271
272

4

Figure 4: Motion comparison of fluids with diﬀerent extension coeﬃcients.
β = 1.005, 1.003, and 1.001 from left to right. Particles are colored according
to their viscosity values (low viscosity: light green, and high viscosity: darkgreen).

al. [52], an elastic force term requires small time steps to simulate elastic materials in a numerically stable manner. To avoid
282 this condition, Zhou et al. [52] relied on an implicit method and283 successfully simulated elastic materials. However, the deforma284 tions of the materials are restricted to some extent because they
285 used a Cauchy tensor, ignoring a second-order term, to build a
286 linear system in an implicit manner. Unlike their method, we
287 still depend on a position-based method to simulate highly de288 formable elastic materials.
289
As the basic concept regarding the handling of elasticity is
290 to correct positions of particles based on constraints, which are
291 the same as viscosity constraints, we also use Eq. (9) as elastic292 ity constraints. For viscous eﬀects, we extend the interparticle
293 distance di j , which is a reference distance when we correct par294 ticle positions, by using Eq. (13). For elasticity, however, we do
295 not extend the interparticle distance di j ; therefore, we can gen296 erate motions of deformable objects with density constraints,
297 which address fluid compression, by restoring the original con298 figurations of particles, as in position-based dynamics [9, 10].
299 We set initial constraints to determine the shape of objects when
300 they are created. If we always correct particle positions using
301 Eqs. (10) and (11) without generating, modifying, and delet302 ing constraints, our method ensures the reproduction of elastic
303 motions of objects.

Figure 5: Two dropped fluid balls mix with each other being cooled by cold
boundary particles in a virtual box on the ground. Each fluid particle has adiﬀerent extension coeﬃcient computed based on particle’s temperature. Particles are colored according to their temperatures (low temperature: blue, midtemperature: green, and high temperature: red).

280
281

309

3.6. Thermal Conduction

The characteristics of materials can change depending on
their temperature. For example, high-temperature materials tend
312 to be less viscous and exhibit liquid-like motions, while low313 temperature materials are highly viscous and can be quasi-rigid
314 or -elastic from fluid through phase transition. To reproduce
315 these phenomena, we incorporate thermal conduction into our
316 framework by computing the temperature of particles on the ba317 sis of thermal diﬀusion equations. Specifically, we improve the318 method proposed by Cleary and Monaghan [53] to handle heat
319 transfer on object surfaces with non-fluid boundary particles:
310
311

1 ∑ 2ci c j m j
T ji ∇2 Wi j
ρi j ci + c j ρ j1 ∑ 2ci ck 1
T ki ∇2 Wik ,+ρi k ci + ck δkTi ← Ti +where T i j = T i − T j and c = sc ρ0 h2 with a parameter sc to
321 control the speed of heat propagation.
322
We then change values in the extension coeﬃcients to vary
323 the viscosity of fluid. Although there are many models for re324 lating temperature and viscosity, we use a simple linear inter325 polation to compute extension coeﬃcient βi with particle tem326 perature for our aim of fast simulation:
304 3.5. Constraint Control AlgorithmβbTb < Ti,305
We summarize control steps for generating, modifying, and
T i −T a
β
+
(β
−
β
)
Ta < Ti ≤ Tb,
β
=

a
b
a
i
306 deleting constraints in Algorithm 1. Note that we use the viscous
T b −T a




β
T
a
i ≤ Ta,307 label for viscous fluid particles and elastic for particles in elastic materials.
327 where βa and βb (βa < βb ) are minimum and maximum val328 ues of βi , respectively, and T a and T b (T a < T b ) are lower
Algorithm 1 Constraint control
329 and upper thresholds of T i , respectively.
This simple model
1: for all fluid particle i do
330 is suﬃcient to generate plausible results, as shown in Figure 5,
2:
for all connected particle j do
331 where a cold ball (blue) and a hot one (red) are dropped onto3:
if at least either i or j is viscous then332 a solid cold ball, and high-temperature particles flow quickly4:
if H ≤ ||pi j || then
333 while low-temperature particles do so slowly. In addition to the5:delete the constraint
334 changes in the extension coeﬃcients, we alter the state of par6:
if not satisfy Eq. (12) then
335 ticles to simulate the phase transition in materials by labeling
7:
modify the constraint
336 the particles as elastic (viscous) if their temperatures are lower
8:
for all neighboring particle j do
337 (higher) than a threshold temperature T c (< T a ).
9:
if at least either i or j is viscous then338
Our method can handle various types of objects, such as10:
if there is no constraint then339 viscous fluids and elastic materials with phase transition and
11:
generate a new constraint
340 changes in physical properties, due to the thermal conduction308
341 in a unified position-based dynamics framework.
320

5

3.7. Procedure of Our Method
Algorithm 2 gives a full-fledged procedure of the steps per344 formed in our method. An object Sn has Mn constraints, de345 noted as C n,1 , . . . , C n,M , which are generated, modified, andn346 deleted using Algorithm 1 (line 11).
We iterate the position
347 correction step until its iteration reaches a specified count. Note
348 that the viscosity constraints are solved with density constraints
349 in parallel to satisfy both types of constraints as much as possi350 ble because the eﬀect of one type of constraint can be ignored
if we solve one side prior to the other.
342

Table 1: Simulation statistics.

343

Algorithm 2 Procedure of our method
1: initialize all variables
2: while animating do
3:for all particle i do
4:
if fluid particle then5:
apply XSPH viscosity
6:
apply external forces vi ← vi + ∆tfext (xi )7:predict position pi ← xi + ∆tvi
8:
update temperature T i ← T i + ∆T i
9:for all particle i do
10:
find neighboring particles at pi
11:control constraints S1 (C1,1 , . . . , C1,M1 ), . . .
12:while correcting particle positions do
13:
for all active non-fluid particle i do
14:
compute δi15:
for all fluid particle i do
16:
compute λi and scorr, i
17:
for all fluid particle i do
18:
compute ∆pdens, i and ∆pvisc, i
19:
for all fluid particle i do
20:
correct position pi ← pi + ∆pdens, i + ∆pvisc, i
21:
for all particle i do
22:
update velocity vi ← (pi − xi )/∆t23:
update position xi ← pi
24:
compute ∆T i
351

352

4. Results

We implemented our algorithm in C++ and parallelized itwith OpenMP 2.0. We used constant time steps in each scene
355 and only six iterations of position correction for fast simulation.
356 In our simulations, parameters were empirically adjusted ac357 cording to the scenes, and non-fluid boundary particles used the358 extension coeﬃcient of their connected particles. All the scenes
359 were executed on a PC with a 4-core Intel Core i7 3.50 GHz
360 CPU and RAM 16.0 GB. We used six threads in total and ob361 tained the accelerated simulation performance by a factor from
362 two to four. Fluid surfaces were extracted with anisotropic ker363 nels [54]. All the results were rendered with a free open-source
364 raytracer POV-Ray 3.7, and oﬀ-line rendering required about 60
365 seconds per frame on average to generate Figure 1. The accom366 panying video shows several results produced with our method
367 and a previous method which uses explicit viscosity integration
368 for comparison. We tabulate the simulation conditions and per369 formances of the results in Table 1.
353
354

Figure #
1
2 (a)
2 (b)
4
5
6 (a)
6 (b)
6 (c)
7
8
9
10
11
12

370

# of particles
fluid/non-fluid
up to 307.2k/111.4k
24.8k/0.0
24.8k/6.8k
27.6k/59.5k
364.1k/126.1k
7.8k/7.7k
7.8k/7.7k
7.8k/7.7k
up to 23.4k/34.0k
up to 92.1k/34.0k
up to 269.4k/111.5k
65.2k/63.7k
up to 33.1k/35.7k
15.9k/75.7k

Time /
frame (s)
47.1
0.5
0.6
1.8
40.5
0.5
82.7
0.5
2.8
6.3
36.2
5.8
1.9
2.4

4.1. Comparison

To prove that our method can use larger time steps thanthose adopted in SPH with explicit viscosity integration, we
373 used a simple scenario in which a viscous fluid ball was dropped
374 onto the ground. As a method for comparison, we used Predictive375 Corrective Incompressible SPH [37], instead of position-based
376 fluids, with Laplacian form of viscosity because the combina377 tion of position-based fluids and Laplacian form of viscosity
378 could not generate plausible viscous fluid behavior. In force379 based methods, using smaller time steps contributes to produc380 ing more accurate simulation results, yet this fact is inapplica381 ble to position-based methods. As noted in the first paper on
382 position-based dynamics proposed by Müller et al. [9], time
383 step size aﬀects simulation results. Specifically, the position384 based dynamics method with smaller time steps generates stiﬀer
385 object motions, while using larger time steps leads to gener386 ating softer and more deformable object motions. Therefore,
387 position-based fluids with extremely small time steps, which
388 are needed for Laplacian form of viscosity, did not produce
389 highly deformable motions of fluid.
Although the compari390 son of the previous method and ours is not completely fair,
391 the viscosity term is a dominant factor for determining time
392 steps in the scene, and this comparison is intended to show
393 that our method can use much larger time steps than the tra394 ditional method to reduce the number of simulation loops each
395 of which includes a time-consuming neighbor search step. Our−4396 method required ∆t = 9.06 × 10
s to produce plausible vis397 cous fluid motions (Figure 6 (a)), while the previous method re−6
398 quired ∆t = 5.41 × 10
s to generate a similar result (Figure 6
399 (b)). As shown in Figure 6 (c), we could not generate a plau400 sible result using the previous method with the same time step
401 used in Figure 6 (a). Our method (Figure 6 (a)) was able to use
402 167.47 times larger time steps than those adopted with the pre403 vious method (Figure 6 (b)); therefore our method could pro404 duce visually plausible viscous fluid motions with much lower
405 computational cost compared to the previous method.
371
372

6

(a) Our method (∆t = 9.06 × 10−4 s)Figure 8: Viscous sheet buckling reproduced with our method.

(b) Laplacian form (∆t = 5.41 × 10−6 s)Figure 9: Several slime lumps with diﬀerent viscosity values successively
dropping onto a solid bunny. Particles are colored according to their viscosity
values (low viscosity: light green, and high viscosity: dark green).

(c) Laplacian form (∆t = 9.06 × 10−4 s)particles used for the positional constraints. In this scene, a viscous lump drops onto a rotating bar, avoids the bar by dividing
425 into two parts, and finally falls onto the ground.
423
424

Figure 6: Comparison of resulting motions.

4.5. Phase Transition
Figure 11 shows the results of a hot viscous fluid poured
428 onto a cold slope. First, the fluid flows on the slope due to low
429 viscosity, while the fluid hardens and finally congeals through430 phase transition as it is cooled by the cold particles of the slope.
431 Figure 12 demonstrates the phase transition in a single object.
432 After an elastic ball runs down a cold slope, the ball bounces
433 a few times, sticks to a hot wall, and starts melting. The com434 bination of viscosity and elasticity constraints can work well
435 through phase transition and generate plausible results.
426
427

Figure 7: Viscous thread coiling reproduced with our method.

406

4.2. Coiling and Buckling

Figures 7 and 8 demonstrate coiling and buckling phenom408 ena, respectively, which highly viscous fluids exhibit. In these409 scenes, a viscous material is dropped onto the ground, while
410 particles are continuously added on top of the material. Our
411 method can generate characteristic and complicated coiling and412 buckling phenomena despite the approximation of viscous fluid
413 dynamics with position-based constraints.
407

4.3. Variable Viscosity
415
Figure 9 provides an example with materials of diﬀerent
416 viscosity values. Several slime lumps with diﬀerent viscosity
417 values are successively dropped onto a solid bunny, mixing with
418 each other. Since each particle has its own parameters, spatial
419 variation in viscosity can be easily achieved.
414

420
421
422

4.4. Boundary-Handling
Figure 10 demonstrates that our boundary-handling scheme
is applicable to a moving solid object with non-fluid boundary436

5. Discussions and Limitations

The configuration of particles is highly relevant to the quality of fluid surfaces due to the dependence of our surface recon439 struction on the particles. Since each particle in position-based
440 fluids has fewer neighboring ones than those in SPH simula441 tions, individual particles significantly aﬀect neighboring ones.
442 Consequently, a particle around a free surface, which originally
443 has fewer neighboring particles, is easily pushed toward the444 surface, causing local irregularities and giving rise to bumpy
445 surfaces. This problem can be alleviated with a sophisticated446 surface reconstruction method [55].
447
Even if viscosity constraints are of no eﬀect, fluid particles
448 on object surfaces are likely to unnaturally adhere to non-fluid
449 ones to satisfy the density constraints by compensating for low
450 particle density.
Although our boundary-handling scheme is451 necessary to specify boundary conditions with constraints, it
452 cannot fully prevent fluid particles from sticking to non-fluid
453 ones. To address this problem, the use of air particles [56] at
454 the significant sacrifice of computational eﬃciency will help us
455 to precisely estimate the particle density on surfaces.
437
438

7

Figure 11: Hot viscous fluid poured onto a cold slope. Lower images are corresponding particle representations, where particles are colored according to theirtemperatures (low temperature: blue, mid-temperature: green, and high temperature: red).

Figure 12: Elastic ball running down a cold slope. Lower images are corresponding particle representations, where particles are colored according to theirtemperatures (low temperature: blue, mid-temperature: green, and high temperature: red).

Our method can encompass a broad range of viscosity: fromalmost inviscid to highly viscous fluids (see Figures 4 and 9).
458 Fortunately, the strength of the viscosity constraints is virtually
459 irrelevant to the required iterations for convergence. We can
460 also use the same time steps regardless of the strength. How461 ever, since time steps can aﬀect the viscosity of fluid, we need
462 to choose appropriate time steps to generate desirable viscous
463 fluid motions. Finding appropriate time steps would be diﬃcult
464 due to the nonlinearity of resolving the constraints and the lack465 of physical grounds, and this point will be investigated further
466 as an important topic for future work.
456
457

associated boundary-handling scheme for position-based fluids appropriately addressed boundary conditions for the con479 straints. Moreover, we incorporated handling of elastic materi480 als into our framework and added thermal conduction to sim481 ulate changes in particle properties and phase transition in the
482 materials. We proved the eﬀectiveness and versatility of our
483 method with several examples.
477
478

484

Acknowledgements

This work has been partly supported by Japan Society for
the Promotion of Science under Grant-in-Aid for Scientific Re487 search (A) No. 26240015. We would like to thank Nobuyuki
467 6. Conclusions
488 Umetani from Disney Research and anonymous reviewers for468
We extended the previous paper [11] and proposed a particle- 489 their valuable suggestions and comments.
485
486

based Lagrangian method for simulating viscous fluids using
position-based dynamics. Fluid volumes are discretized by par471 ticles that interact with other neighboring ones under position472 based constraints, approximating the dynamics of viscous flu473 ids.
Hence, our method allowed for the use of larger time
474 steps than those adopted in SPH methods with explicit viscosity
475 integration, while enabling us to reproduce plausible motions
476 of fluids with diﬀerent properties in a unified framework. An
469
470

490

491
492
493
494
495
496
497

8

References[1] Carlson M, Mucha PJ, Van Horn III RB, Turk G. Melting and flowing. In:
Proceedings of the 2002 ACM SIGGRAPH/Eurographics Symposium on
Computer Animation. 2002, p. 167–74.
[2] Rasmussen N, Enright D, Nguyen D, Marino S, Sumner N, Geiger W,
et al. Directable photorealistic liquids. In: Proceedings of the 2004 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation. 2004, p.
193–202.

Astronomy and Astrophysics 1992;30:543–74.
[20] Stora D, olivier Agliati P, paule Cani M, Neyret F, dominique Gascuel J.
545
Animating lava flows. In: Graphics Interface. 1999, p. 203–10.
546 [21] Solenthaler B, Schläfli J, Pajarola R. A unified particle model for fluid
547
solid interactions: Research articles. Computer Animation and Virtual
548
Worlds 2007;18(1):69–82.
549 [22] Becker M, Ihmsen M, Teschner M. Corotated SPH for deformable solids.
550
In: Proceedings of Eurographics Workshop on Natural Phenomena. 2009,551
p. 27–34.
552 [23] Monaghan JJ. On the problem of penetration in particle methods. Journal553
of Computational Physics 1989;82(1):1–15.
554 [24] Mao H, Yang YH. Particle-based non-newtonian fluid animation with555
heating eﬀects. Tech. Rep.; University of Alberta; 2006.
556 [25] Chang Y, Bao K, Liu Y, Zhu J, Wu E. A particle-based method for vis557
coelastic fluids animation. In: Proceedings of the 16th ACM Symposium558
on Virtual Reality Software and Technology. 2009, p. 111–7.
559 [26] Rafiee A, Manzari M, Hosseini M. An incompressible SPH method for
560
simulation of unsteady viscoelastic free-surface flows. International Jour561
nal of Non-Linear Mechanics 2007;42(10):1210 –23.
562 [27] Gerszewski D, Bhattacharya H, Bargteil AW.
A point-based method563
for animating elastoplastic solids. In: Proceedings of the 2009 ACM
564
SIGGRAPH/Eurographics Symposium on Computer Animation. 2009, p.
565
133–8.
566 [28] Takamatsu K, Kanai T. A fast and practical method for animating particle567
based viscoelastic fluids. The International Journal of Virtual Reality568
2011;10:29–35.
569 [29] Dagenais F, Gagnon J, Paquette E. A prediction-correction approach for
570
stable SPH fluid simulation from liquid to rigid. In: Proceedings of the571
Computer Graphics International 2012. 2012,.
572 [30] Miller G, Pearce A. Globular dynamics: A connected particle system for573
animating viscous fluids. Computers and Graphics 1989;13(3):305–9.
574 [31] Terzopoulos D, Platt J, Fleischer K. Heating and melting deformable
575
models. Journal of Visualization and Computer Animation 1991;2:68–
576
73.
577 [32] Steele K, Cline D, Egbert PK, Dinerstein J. Modeling and rendering vis578
cous liquids. Computer Animation and Virtual Worlds 2004;15(3-4):183–579
92.
580 [33] Tamura N, Nakaguchi T, Tsumura N, Miyake Y. Spring-bead anima581
tion of viscoelastic materials. IEEE Computer Graphics and Applications
582
2007;27(6):87–93.
583 [34] Clavet S, Beaudoin P, Poulin P. Particle-based viscoelastic fluid sim584
ulation. In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics
585
Symposium on Computer Animation. 2005, p. 219–28.
586 [35] Ihmsen M, Orthmann J, Solenthaler B, Kolb A, Teschner M. SPH fluids in
587
computer graphics. In: EUROGRAPHICS 2014 State of the Art Reports.
588
2014, p. 21–42.
589 [36] Becker M, Teschner M. Weakly compressible SPH for free surface flows.
590
In: Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium591
on Computer Animation. 2007, p. 209–17.
592 [37] Solenthaler B, Pajarola R.
Predictive-corrective incompressible SPH.
593
ACM Transactions on Graphics 2009;28(3):1–6.
594 [38] Ihmsen M, Cornelis J, Solenthaler B, Horvath C, Teschner M. Implicit
595
incompressible SPH. IEEE Transactions on Visualization and Computer
596
Graphics 2014;20(3):426–35.
597 [39] Becker M, Tessendorf H, Teschner M. Direct forcing for Lagrangian598
rigid-fluid coupling. IEEE Transactions on Visualization and Computer
599
Graphics 2009;15(3):493–503.
600 [40] Ihmsen M, Akinci N, Gissler M, Teschner M. Boundary handling and601
adaptive time-stepping for PCISPH. In: Proceedings of 10th Workshop
602
on Virtual Reality Interaction and Physical Simulation. 2010, p. 79–88.
603 [41] Akinci N, Ihmsen M, Akinci G, Solenthaler B, Teschner M.
Versa604
tile rigid-fluid coupling for incompressible SPH. ACM Transactions on
605
Graphics 2012;31(4):1–8.
606 [42] Akinci N, Cornelis J, Akinci G, Teschner M. Coupling elastic solids with607
smoothed particle hydrodynamics fluids. Computer Animation and Vir608
tual Worlds 2013;24(3-4):195–203.
609 [43] Müller M, Solenthaler B, Keiser R, Gross M. Particle-based fluid-fluid
610
interaction. In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics
611
Symposium on Computer Animation. 2005, p. 237–44.
612 [44] Solenthaler B, Pajarola R. Density contrast SPH interfaces. In: Proceed613
ings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Com543
544

Figure 10: Viscous lump dropping onto a rotating bar. Lower images arecorresponding particle representations.

[3] Batty C, Bridson R. Accurate viscous free surfaces for buckling, coiling, and rotating liquids. In: Proceedings of the 2008 ACM SIG500
GRAPH/Eurographics Symposium on Computer Animation. 2008, p.
501
219–28.
502 [4] Batty C, Houston B. A simple finite volume method for adaptive viscous
503
liquids. In: Proceedings of the 2011 ACM SIGGRAPH/Eurographics
504
Symposium on Computer Animation. 2011, p. 111–8.
505 [5] Paiva A, Petronetto F, Lewiner T, Tavares G. Particle-based viscoplastic
506
fluid/solid simulation. Computer-Aided Design 2009;41(4):306–14.
507 [6] Müller M, Charypar D, Gross M.
Particle-based fluid simulation for
508
interactive applications. In: Proceedings of the 2003 ACM SIG509
GRAPH/Eurographics Symposium on Computer Animation. 2003, p.
510
154–9.
511 [7] Charney JG, Fjortoft R, Von Neumann J. Numerical integration of the512barotropic vorticity equation. Tellus 1950;2(4):237–54.
513 [8] Macklin M, Müller M. Position based fluids. ACM Transactions on
514
Graphics 2013;32(4):1–5.
515 [9] Müller M, Heidelberger B, Hennix M, Ratcliﬀ J. Position based dy516
namics. Journal of Visual Communication and Image Representation
517
2007;18(2):109–18.
518 [10] Bender J, Müller M, Otaduy MA, Teschner M. Position-based methods
519
for the simulation of solid objects in computer graphics. In: EURO520
GRAPHICS 2013 State of the Art Reports. 2013, p. 1–22.
521 [11] Takahashi T, Fujishiro I.
Accelerated viscous fluid simulation using522
position-based constraints. In: Proceedings of CAD/Graphics 2013. 2013,523
p. 260–7.
524 [12] Bargteil AW, Wojtan C, Hodgins JK, Turk G. A finite element method
525
for animating large viscoplastic flow. ACM Transactions on Graphics
526
2007;26(3).
527 [13] Wojtan C, Turk G. Fast viscoelastic behavior with thin features. ACM
528
Transactions on Graphics 2008;27(3):1–8.
529 [14] Wicke M, Ritchie D, Klingner BM, Burke S, Shewchuk JR, O’Brien JF.
530
Dynamic local remeshing for elastoplastic simulation. ACM Transactions
531
on Graphics 2010;29:1–11.
532 [15] Clausen P, Wicke M, Shewchuk JR, O’Brien JF. Simulating liquids and533
solid-liquid interactions with Lagrangian meshes. ACM Transactions on
534
Graphics 2013;32(2):1–15.
535 [16] Bergou M, Audoly B, Vouga E, Wardetzky M, Grinspun E. Discrete vis536
cous threads. ACM Transactions on Graphics 2010;29(4):1–10.
537 [17] Batty C, Uribe A, Audoly B, Grinspun E. Discrete viscous sheets. ACM
538
Transactions on Graphics 2012;31(4):1–7.
539 [18] Desbrun M, Gascuel MP. Smoothed particles: A new paradigm for an540
imating highly deformable bodies. In: Proceedings of the Eurographics
541
Workshop on Computer Animation and Simulation. 1996, p. 61–76.
542 [19] Monaghan JJ. Smoothed particle hydrodynamics. Annual Review of
498
499

9

puter Animation. 2008, p. 211–8.
[45] Orthmann J, Hochstetter H, Bader J, Bayraktar S, Kolb A. Consistent616
surface model for SPH-based fluid transport. In: Proceedings of the 2013
617
ACM SIGGRAPH/Eurographics Symposium on Computer Animation.
618
2013, p. 95–103.
619 [46] Ihmsen M, Akinci N, Becker M, Teschner M. A parallel SPH implemen620
tation on multi-core CPUs. Computer Graphics Forum 2011;30(1):99–621
112.
622 [47] Harada T, Koshizuka S, Kawaguchi Y. Smoothed particle hydrodynamics
623
on GPUs. In: Proceedings of Computer Graphics International. 2007, p.
624
63–70.
625 [48] Goswami P, Schlegel P, Solenthaler B, Pajarola R. Interactive SPH sim626
ulation and rendering on the GPU. In: Proceedings of the 2010 ACM
627
SIGGRAPH/Eurographics Symposium on Computer Animation. 2010, p.
628
55–64.
629 [49] Adams B, Pauly M, Keiser R, Guibas LJ. Adaptively sampled particle
630
fluids. ACM Transactions on Graphics 2007;26(3).
631 [50] Solenthaler B, Gross M. Two-scale particle simulation. ACM Transac632
tions on Graphics 2011;30(4):1–8.
633 [51] Orthmann J, Kolb A. Temporal blending for adaptive SPH. Computer
634
Graphics Forum 2012;31(8):2436–49.
635 [52] Zhou Y, Lun Z, Kalogerakis E, Wang R. Implicit integration for particle636
based simulation of elasto-plastic solids. Computer Graphics Forum637
2013;32(7):215–23.
638 [53] Cleary PW, Monaghan JJ. Conduction modelling using smoothed particle
639
hydrodynamics. Journal of Computational Physics 1999;148(1):227–64.
640 [54] Yu J, Turk G.
Reconstructing surfaces of particle-based fluids using641
anisotropic kernels. ACM Transactions on Graphics 2013;32(1):1–12.
642 [55] Ando R, Thürey N, Wojtan C. Highly adaptive liquid simulations on
643
tetrahedral meshes. ACM Transactions on Graphics 2013;32(4):1–10.
644 [56] Schechter H, Bridson R. Ghost SPH for animating water. ACM Transac645
tions on Graphics 2012;31(4):1–8.
614
615

10