Noname manuscript No.
(will be inserted by the editor)

Volume Preserving Viscoelastic Fluids with Large
Deformations Using Position-based Velocity Corrections
Tetsuya Takahashi · Yoshinori Dobashi · Issei Fujishiro · Tomoyuki
Nishita

Received: date / Accepted: date

Abstract We propose a particle-based hybrid method
for simulating volume preserving viscoelastic fluids withlarge deformations. Our method combines Smoothed Particle Hydrodynamics (SPH) and Position-based Dynamics (PBD) to approximate the dynamics of viscoelastic fluids. While preserving their volumes using SPH,
we exploit an idea of PBD and correct particle velocities for viscoelastic eﬀects not to negatively aﬀect volume preservation of materials. To correct particle velocities and simulate viscoelastic fluids, we use connections
between particles which are adaptively generated anddeleted based on the positional relations of the particles. Additionally, we weaken the eﬀect of velocity corrections to address plastic deformations of materials.
For one-way and two-way fluid-solid coupling, we incorporate solid boundary particles into our algorithm.
Several examples demonstrate that our hybrid method
can suﬃciently preserve fluid volumes and robustly andplausibly generate a variety of viscoelastic behaviors,such as splitting and merging, large deformations, andBarus eﬀect.
Keywords Fluid simulation · Viscoelasticity ·Deformation · Volume preservation · Velocity correction
Tetsuya Takahashi
UEI Research / Keio University
E-mail: tetsuya.takahashi@uei.co.jp
Yoshinori Dobashi
Hokkaido University / UEI ResearchE-mail: doba@ime.ist.hokudai.ac.jp
Issei Fujishiro
Keio University
E-mail: fuji@ics.keio.ac.jp
Tomoyuki Nishita
UEI Research / Hiroshima Shudo University
E-mail: tomoyuki.nishita@uei.co.jp

1 Introduction
Particle-based methods have been popular techniques
to simulate fluids, rigid bodies, deformable objects, andtheir interactions. These methods use a simple conceptof particle interactions, which can be easily extended tosimulate various types of objects, and were adopted fora versatile simulation framework [16]. For conceptual
simplicity and versatility, we herein focus on particlebased methods.
Since viscoelastic materials play an important role
in representing common materials (e.g., egg white, gels,
toothpaste, and slime) and producing visually attractive characters and eﬀects, simulating viscoelastic materials has been required for movies and video games. In
such entertainments, exaggerated representations whichreal materials do not exhibit, e.g., very large deformations, are frequently required to make characters andeﬀects appear interesting and impressive.
In computer graphics, particle-based methods for
simulating viscoelastic fluids have been proposed. However, previous particle-based, Smoothed Particle Hydrodynamics (SPH) methods which depend on physicallybased viscoelastic models [17, 6] cannot simulate highly
deformable objects because, at a particle position, physical contributions from farther particles are likely tobe smaller owing to the decreasing property of finite
support kernels, and consequently, viscoelastic materials fail to restore their original shapes when they undergo significant deformations. This problem can be addressed by combining SPH and a geometric method[28, 7]. However, these hybrid methods have anotherproblem that they suﬀer from volume loss of materials,
which leads to unphysical fluid behaviors, because thegeometric methods adopted by Takamatsu and Kanai

2

Tetsuya Takahashi et al.

Fig. 1 Viscoelastic fluids simulated with our hybrid method. (Left) Merging spheres with diﬀerent viscoelasticity values
colliding with a static solid dragon. (Middle) A viscoelastic ball colliding with a solid cube. (Right) Barus eﬀect of a viscoelastic
material (see §4.4).

[28] and Clavet et al. [7] make an incompressible fluid
solver fail to preserve fluid volumes.
To solve the two problems above simultaneously, we
propose a particle-based hybrid method for simulating
viscoelastic fluids, which can preserve fluid volumes andenable large deformations of the fluids. Our method
uses a combination of SPH and Position-based Dynamics (PBD) [20, 4]. While using SPH to preserve fluid
volumes, we exploit PBD for velocity corrections which
achieve viscoelastic eﬀects involving large deformations
without negatively influencing the volume preservation.
To correct particle velocities, we use a set of connections
between particles, which are adaptively updated basedon the positional relations of the particles. We incorporate solid boundary particles into our algorithm tosimulate adhesion of viscoelastic fluids to solid objects.
Our key contribution lies in velocity corrections whichallow for combining a geometric method with SPH to
simulate viscoelastic behavior with large deformations
while preserving fluid volumes, since the SPH-based
methods [17, 6] cannot handle large deformations andthe previously proposed hybrid methods [28,7] cannot
preserve the volume of viscoelastic fluids. Fig. 1 illustrates characteristic behaviors of viscoelastic materials
simulated with our hybrid method.
This paper is an extended version of our previous
paper [26]. In this paper, we focus on the positionbased velocity correction scheme. As new additional
materials, we include a comparison of our method with
the SPH-based method [17] and propose a boundaryhandling method for fluid-solid coupling.

2 Related Work
Viscoelastic fluids Miller and Pearce [18] and Terzopoulos et al. [30] proposed a spring-based model thatcomputes repulsion and attraction forces among particles to simulate viscoelastic materials. Their method

was adopted by Steele et al. [25] and Tamura et al. [29].
Clavet et al. [7] combined this spring-based methodwith SPH to simulate materials with elasticity, plasticity, and viscosity, adopting a prediction-relaxation
scheme. A similar spring-based method was also proposed by Takahashi et al. [27], who used PBD to simulate fluids with viscosity and elasticity in a unified
framework. Takamatsu and Kanai [28] combined ShapeMatching (SM), which was originally proposed by Müller
et al. [21] to simulate deformable objects, with SPH to
fast and robustly simulate viscoelastic fluids.
Müller et al. [22] proposed an elasticity term whichuses Moving Least Square (MLS) to simulate elastoplastic objects. Solenthaler et al. [24] adopted the formulation of this elasticity term and computed it using SPH instead of MLS to allow for robustly simulating fluid with various properties under the condition ofcollinear or coplanar particle distributions. The methodof Solenthaler et al. [24] was extended to handle rotational motions of elastic materials [2]. Mao and Yang[17] introduced a viscoelastic force term, called nonlinear corotational Maxwell model, into the Navier-Stokes
equations to simulate viscoelastic fluids. A method similar to the method of Mao and Yang [17] was proposedby Chang et al. [6].
Gerszewski et al. [9] proposed a new formulation
for simulating elastoplastic materials, which uses aﬃnetransformations to compute the gradient of deformations. This formulation was solved by Zhou et al. [31]in an implicit manner to eﬃciently and robustly perform simulations with larger time steps. The method ofGerszewski et al. [9] was extended by Jones et al. [14].
Jones et al. [14] used MLS to robustly simulate elastoplastic materials with spatially and temporally varying
masses.
Volume preservation In particle-based simulations,volume preservation has been achieved by minimizing
density fluctuations from the fluid rest density by en-

Volume Preserving Viscoelastic Fluids with Large Deformations Using Position-based Velocity Corrections

forcing fluid incompressibility [13] since meshes which
are often used for volume preservation of deformableobjects are inappropriate because of frequent topology changes. After Müller et al. [19] presented the basic SPH method, which suﬀers from unacceptable density fluctuations and volume loss because of the dependence on an equation of state (EOS), Becker andTeschner [3] proposed Weakly Compressible SPH andalleviated large density fluctuations by replacing the
EOF in [19] with a stiﬀer EOS called Tait equation,requiring exceedingly smaller time steps for stable simulation. Solenthaler et al. [23] proposed a predictivecorrective scheme, called Predictive-Corrective Incompressible SPH (PCISPH) to use larger time steps whilesatisfying fluid incompressibility. A similar predictivecorrective scheme called local Poisson SPH was proposed by He et al. [10]. Bodin et al. [5] proposed asystem of velocity constraints to improve the accuracyof PCISPH [23], while achieving incompressible flows.
To further accelerate enforcing fluid incompressibility
by allowing us to use larger time steps than those usedin PCISPH [23], Macklin and Müller [15] adapted PBD
for fluid simulation, and Ihmsen et al. [12] proposed
Implicit Incompressible SPH (IISPH), which computes
pressure values using semi-implicit integration. A methodthat enforces fluid incompressibility using semi-implicit
integration, called Incompressible SPH (ISPH) was also
proposed by Cummins and Rudman [8].

ticle i has position xi , velocity vi , and mass mi ), andpressure pi is iteratively corrected to obtain pressure
force Fpi which achieves small deviations of density ρito the rest density ρ0 less than a certain criterion, e.g.,
0.01%. We briefly explain the flow of computations forthe solvers. For detailed simulation procedures and discussions, refer to the papers of particle-based incompressible fluid solvers [8, 23, 12].
We first compute intermediate velocity vi∗ for particle i byvi∗ = vit + ∆t

Fotheri,mi(1)

where t denotes time, ∆t time step, and Fother all forces
except pressure and viscoelastic ones. Next, we use an
incompressible fluid solver to obtain pressure pi , which
is clamped as pi = 0 if pi < 0 not to cause attraction
forces between particles, through iterative corrections
with intermediate velocity vi∗ and compute pressure
force Fpi by
)
(
∑
pj
pi
p
Fi = −mi+ 2 ∇Wij ,
(2)
mj
ρ2i
ρj
j
where j denotes a neighboring fluid particle and Wij akernel with a kernel radius h. Then, we compute velocity vi for incompressible flows:
vit+1 = vi∗ + ∆t

3 Proposed Method

3

Fp
.
mi

(3)

Finally, we update particle position xi :The key of our method is to separately achieve volume
preservation and large deformations of viscoelastic materials. We rely on a particle-based incompressible fluid
solver to preserve fluid volumes while exploiting an idea
of PBD [20, 4] to correct particle velocities for viscoelastic eﬀects without negatively aﬀecting convergence ofdensity fluctuations in the fluid solver.
We briefly describe computational process and important features of particle-based incompressible fluid
solvers in § 3.1 and explain our velocity correction scheme,
which enables large deformations of viscoelastic materials while preserving their volumes in § 3.2. In § 3.3,
we summarize our algorithm, giving a list of simulation
procedures and implementation details.

3.1 Particle-based Incompressible Fluid Solvers
To preserve fluid volumes, we use a particle-based incompressible fluid solver, such as ISPH [8], PCISPH[23], and IISPH [12]. In these SPH-based incompressible fluid solvers, fluids are discretized by particles (par-xt+1= xti + ∆tvi .
i(4)

There are two important requirements to note; first,
particle velocity vi needs to be updated with pressure
force Fpi with Eq. (3) to integrate particle position xi
with Eq. (4), achieving incompressible flows. However,
the method of Takamatsu and Kanai [28] does not satisfy this requirement because their method mixes velocities derived from SM, which does not ensure incompressible flows, and SPH, and thus resulting velocities
cannot achieve incompressible flows. Second, particle
position xi must be fixed after neighbor search steps
until pressure force Fpi is computed with Eq. (2). As
explained in [23, 12], neglecting changes of particle positions causes erroneous sets of neighboring particles andtheir interparticle distances and leads to inaccurate estimates of physical quantities, which can make the fluid
solver fail to converge or delay the convergence of theiterations. Since the method of Clavet et al. [7] uses
a prediction-relaxation scheme, which change particle
positions, this requirement is not satisfied.

4

Tetsuya Takahashi et al.

3.2 Velocity Correction Scheme

Fig. 2 A ball with diﬀerent viscoelasticity values dropped
onto the ground. Particles are color coded based on theirvelocity correction coeﬃcients c (low: cyan and high: darkcyan).

Influence for viscoelastic effects

Satisfying the two requirements mentioned above, we
correct particle velocities to achieve viscoelastic eﬀects
with velocity correction vector ∆vi (see § 3.2.1). We
obtain intermediate velocity vi∗ by modifying Eq. (1):Fothervi∗ = vit + ∆t i
+ ∆vi .
(5)mi
In [26], although three velocity correction schemes(shape matching, spring-based, and position-based) were
proposed to compute velocity correction vector ∆v, weherein use the position-based scheme because of thefollowing advantages over the other schemes; first, theposition-based scheme requires fewer neighbor particles
than the shape matching scheme and uses pairwise connections which can simplify boundary conditions withother objects unlike the shape matching scheme whichneeds to construct particle clusters including farther
particles. Second, the position-based scheme inherits
properties of PBD [20, 4] and thus is robust and easy
to tune parameters as compared to the spring-based
scheme.

SPH-based method
Our methodrijh||xij||

3.2.1 Position-based Velocity Correction
To compute velocity correction vector ∆v, we construct
a set of pairwise connections for neighboring particles
when an object is created, with their initial interparticle
distance rij , which is initialized as rij = ||xij ||, where
xij = xi − xj . Hereafter, we call particles connectedwith other particles as connected particles. We correct
velocities of connected particles only when their interparticle distance ||xij || is larger than their initial interparticle distance rij (fluid expansion), namely rij <||xij ||, because fluid compression, which generally occurs when ||xij || < rij , is resolved to preserve fluid volumes with our fluid solver. To compute ∆v for thoseconnected particles such that rij < ||xij ||, we define a
distance function Dij as
Dij = max(||xij || − rij , 0).
Then, we compute ∆vi with a velocity correction coeﬃcient ci (0 ≤ li ≤ ci ), where li is the lower limit ofci :Fig. 3 Comparison of the influence for viscoelastic eﬀects
between an SPH-based method and our method.

Our velocity correction vector is basically obtained by
dividing ∆xi by ∆t. By using Dij , the velocity correction vector becomes zero when ||xij || < rij . We additionally introduce velocity correction coeﬃcients ci andcj to control the stiﬀness of viscoelastic materials; lower(higher) ci leads to softer (stiﬀer) materials (see Fig. 2).
Our velocity corrections enable large deformations
that cannot be generated with SPH-based methods,e.g., [17, 6], because the influence of the SPH-based methods for viscoelastic eﬀects decreases as interparticle distances get larger and finally becomes zero owing to finite support kernels when h < ||xij || while the influenceof our method increases when rij < ||xij || (see Fig. 3).

Sf

i
1 ∑
mj
xij
ci + cj∆vi = −Dij,∆t j
2 mi + mj||xij ||

(6)

where Sif denotes the number of connected fluid particles to particle i. To derive Eq. (6), we adopted an
idea of position correction used in PBD [20,4], where
position correction vector is defined asf∆xi = −

Si
∑
j

xij
mj
(||xij || − rij ).
mi + mj||xij ||

3.2.2 Coeﬃcient Relaxation
When large stress is applied to viscoelastic materials,
they cannot completely return to their original shapes(this is known as plastic deformation). To simulate thisbehavior, we weaken the eﬀect of ci if li < ci , taking
into account the fact that connected particles do notrestore their positional relation when their connection
is extended larger than a certain distance. Specifically,Volume Preserving Viscoelastic Fluids with Large Deformations Using Position-based Velocity Corrections

5

3.2.4 Boundary-handling for Solids
We simulate interactions of viscoelastic fluids and solids
by incorporating solid boundary particles into our algorithm. Connections between fluid and solid particles
are updated following Algorithm 1. We compute ∆vi
by extending Eq. (6) with solid particles:
Sf

Fig. 4 (Top) Two merging viscoelastic balls. (Bottom) A
viscoelastic ball thrown toward a column, splitting into threelumps.

ici + cj
1 ∑
mj
xij
Dij
∆vi = −
∆t j
2 mi + mj||xij ||

Ss

i
1 ∑
mk
xik
ci + ck−Dik,∆t
2 mi + mk||xik ||

(8)k

we update ci byf= max(cti − ∆t
ct+1
i

Si
∑
wi + wj

{
σij =

j

2
γ +γ

σij , li ),1 if i 2 j <0 otherwise

Dij
h

(7)

,

where k denotes a neighboring solid particle and Sis thenumber of connected solid particles to particle i. Since
force-based rigid body simulators are often adopted, we
compute viscoelastic force Fvi to exert the antisymmetric eﬀect of fluid particles to solid particle i, preserving
their momentum:fmi
Fvi = − 2∆twhere wi (0 ≤ wi ) controls the weakening speed of ci ,and γi (0 ≤ γi ) is a yield criterion. σij works so thatthe weakening happens only when normalized extension
γ +γ
Dij /h is larger than i 2 j .

Si∑
ci + cj

2

j

mj
xij
Dij
.
mi + mj||xij ||

(9)Computing velocity correction vector ∆vi for oneway solid-to-fluid coupling is straightforward. Assuming mk = ∞ in Eq. (8) for a solid particle k, which isfixed or moves along a determined path, we obtainSf

3.2.3 Connection Control
We adaptively generate and delete particle connections
depending on the positional relations of particles to allow viscoelastic fluids to merge with others and split
into several lumps. For merging, we generate a newconnection between two particles with αi (0 < αi ) if
α +α||xij ||/h < i 2 j , and if there is no connection between
them. For splitting, on the other hand, we remove partiβ +β
cle connections with βi (0 < βi ) if i 2 j < ||xij ||/h. We
show our connection control algorithm in Algorithm 1for clarity. Merging and splitting of viscoelastic materials are demonstrated in Fig. 4.

ici + cj
1 ∑mjxij∆vi = −Dij∆t j
2 mi + mj||xij ||

Ssi
1 ∑ci + ck
xik
−
Dik
.
∆t
2
||xik ||

(10)k

Extending Eq. (7) with solid particles, we weaken
the eﬀect of ci asf

ct+1= max(cti − ∆t(
i

Si
∑
j

s

wij σij +Si∑

wik σik ), li ),

(11)

k

w +w

where wij = i 2 j .
Fig. 5 illustrates the eﬀects of the fluid-solid coupling with an example of a cube thrown toward a wall.

Algorithm 1 Connection control algorithm
1: for all fluid particle i do2:
for all neighboring particle j do
α +α3:
if ||xij ||/h < i 2 j then4:
if there is no connection then5:generate a new connection between i and j
6:
for all connected particle j do
β +β7:
if i 2 j < ||xij ||/h then8:
delete the connection between i and j

3.3 Simulation Algorithm and Implementation
We give an outline of our procedures in Algorithm 2.
In our implementation, we use IISPH [12] as an incompressible fluid solver and set a convergence criterion as 0.01%. For fluid-solid coupling of pressure andviscosity, we employ the method of Akinci et al. [1].
We use finite support kernels as proposed by Müller et6Tetsuya Takahashi et al.
Table 1 Simulation conditions and performance. N f and N sdenote the number of fluid and solid particles, respectively.
ttotal is total simulation time per frame.

Fig. 5 A cube thrown toward a wall without (left) / with(right) fluid-solid coupling.

al. [19]. To accelerate neighbor search steps, we used avariant of the z-index sort algorithm presented in [11].
There are six adjustable parameters c, l, w, γ, α, and βthat we introduced to cover various phenomena of viscoelastic fluids, and we basically use values similar toc = 0.001, l = 0.0001, w = 0.001, γ = 1.01, α = 1.0, and
β = 2.0, adjusting them according to scenarios.

Figure #2 (left)
2 (middle)
2 (right)
4 (top)4 (bottom)
5 (left)
5 (right)
6 (top)6 (bottom)
7 (middle)
7 (right)
9 (middle)
9 (right)
1112
13

N f /N s
24.5k/38.6k
24.5k/38.6k
24.5k/38.6k
45.2k/60.3k
22.6k/67.4k
13.8k/27.0k
13.8k/27.0k
7.1k/21.8k
7.1k/21.8k
8.2k/19.2k
8.2k/19.2k
14.0k/27.5k
14.0k/27.5k
up to 195.1k/168.1k
33.4k/65.7k
73.5k/108.8k

ttotal (s)
3.02
3.83
4.06
9.57
5.36
2.40
2.43
0.39
0.64
5.75
2.47
3.35
3.13
125.2
6.82
29.01

Algorithm 2 Simulation algorithm
1: for all particle i do
2:find neighboring particles
3: for all particle i do
4:
update particle connections following Algorithm 1
5:
if particle i is a fluid particle then
6:obtain vi∗ with Eqs. (5), (8), and (10)7:relax ci with Eq. (11)8: compute particle pressure p using IISPH
9: for all fluid particle i do
10:compute Fp
i with Eq. (2)update solids with Eq. (9)11: for all fluid particle i do
12:
integrate vi with Eq. (3)13:integrate xi with Eq. (4)4 Results
We implemented our method in C++ and parallelized itusing OpenMP 2.0. All the simulations were performed
on a PC with a 4-core Intel Core i7 3.40 GHz CPU
and RAM 16.0 GB. We used a physically-based renderer Mitsuba to render all the scenarios. Our velocity
corrections generally occupy only 1.3% of whole simulation time. Simulation conditions and performance arelisted in Table 1.

4.1 Large Deformation
We compare results generated with our method andthe SPH-based method using a non-linear corotational
Maxwell model proposed by Mao and Yang [17] to showthat our hybrid method can handle large deformationsof a viscoelastic material, which cannot be simulatedwith the SPH-based method [17]. Though the methodof Mao and Yang [17] resampled particles to solve anissue of particle disorder, we did not do that in our
experiment to show that our method can generate large
deformations without the need of special treatments.
Fig. 6 illustrates a viscoelastic sphere dropped onto
the ground. The viscoelastic ball simulated with themethod of Mao and Yang [17] exhibits slight elastic
deformations. After the degree of the deformations exceeds a certain threshold, the ball fails to restore their
original shape and breaks into many lumps which consist of fairly close particles because particles close toothers exert strong forces to preserve their positional
relations (see Fig. 3). The method of Mao and Yang[17] needs to use strong viscoelastic forces to make theviscoelastic ball exhibit slight elastic deformations at
the expense of numerical stability. On the other hand,
our method can handle large deformations, restoring its
original shape without resampling particles nor introducing numerical instability.

4.2 Volume Preservation
We compare our method and the existing hybrid methods [28, 7], using a scenario of a viscoelastic ball compressed by a moving board along a vertical path.
Comparison to [28] First, we illustrate the diﬀerence between our method and the method of Takamatsu
and Kanai [28], who combined SPH and SM. Since thebasic SPH [19] adopted by Takamatsu and Kanai [28] astheir underlying fluid solver cannot preserve fluid vol-

Volume V (m3)Volume Preserving Viscoelastic Fluids with Large Deformations Using Position-based Velocity Corrections
Takamatsu and Kanai
Our method
Reference

Ϭ͘Ϭϵϱ
Ϭ͘Ϭϴϱ
Ϭ͘Ϭϳϱ
Ϭ͘Ϭϲϱ
Ϭ͘Ϭϱϱ
Ϭ͘Ϭϰϱ
Ϭ͘Ϭϯϱ
Ϭ͘ϬϮϱ
Ϭ͘Ϭϭϱ
Ϭ

Fig. 6 Comparison for large deformation. (Top) the methodof Mao and Yang [17]. (Bottom) our method.

Fig. 7 Comparison for volume preservation. (Left) initialstate. (Middle) the method of Takamatsu and Kanai [28].
(Right) our method.

umes, we use IISPH [12] instead of the basic SPH [19]and consider the method of [28] as a hybrid of IISPHand SM. In this comparison, we use 100 pressure iterations in IISPH for both methods because the methodof [28] is likely to fail to preserve fluid volumes even if
more pressure iterations are applied.
Fig. 7 shows the comparison, where the left is theinitial state, and the middle (right) is the result of themethod [28] (our method), and Fig. 8 illustrates thevolume V (which is estimated by summingup∑ the vol∑ume of all fluid particles as V = i Vi = i mi /ρi )
occupied by the ball for the method of Takamatsu andKanai [28], our method, and reference (initial volume).
The volume of the ball on the middle decreases and increases depending on the moving board. By contrast,
the ball on the right preserves its volume, which isfairly close to the reference. Additionally, when being
compressed, the ball on the right exhibits horizontally
spreading motions, which cannot be generated with themethod of [28].

7

ϱϬ

ϭϬϬ
FramesϭϱϬϮϬϬFig. 8 Profiles of the volume V occupied by a ball for Fig. 7

Fig. 9 Comparison for volume preservation. (Left) initialstate. (Middle) the method of Clavet et al. [7] (Right) our
method.

Comparison to [7] Second, we compare our method
and the method of Clavet et al. [7], who combined
SPH with spring systems using a prediction-relaxation
scheme. Similar to the above previous method, since
Clavet et al. [7] used the basic SPH [19], we use IISPH
as their fluid solver for the same reason.

Fig. 9 demonstrates the comparison, where the left
is the initial state, and the middle (right) is the resultof the method [7] (our method). Since errors introducedby position changes during the iteration are relatively
small and it is diﬃcult to clarify visual or volume diﬀerences between the previous method [7] and our method,
we compare these methods in terms of the number of
iterations required to satisfy the convergence criterionof 0.01%, as shown in Fig. 10, where profiles of iteration
numbers for Clavet. et al [7] nspr and our method npos ,and their diﬀerence nspr −npos are illustrated. Although
our method and the method of Clavet et al. [7] can genMoreover, our method is advantageous over the method erate similar results, the method of Clavet et al. [7] isof [28] in terms of performance (see Table 1), and there
likely to require more iterations to enforce fluid incomare two main factors; first, the method of Takamatsupressibility than our method. This is noticeable espeand Kanai [28] needs more neighboring particles to stacially when the ball is compressed and significantly debilize simulations than our method, and processing these formed by the board because large deformations cause
particles is an additional cost. Second, their method [28]strong attraction forces and introduce larger positional
uses Singular Value Decomposition, which is time conerrors (see “Diﬀerence” in Fig. 10). The extra iterations
suming, to obtain rotational matrices in SM while ourof the method of Clavet et al. [7] lead to slightly higher
velocity corrections need simple computations only.
computational cost than our method (see Table 1).

8

Tetsuya Takahashi et al.

Clavet et al.
Our method
Difference

Iteration number

ϭϯϱ
ϭϭϱ
ϵϱ
ϳϱ

Fig. 13 Barus eﬀect of a viscoelastic material.

ϱϱ
ϯϱ
ϭϱ
Ͳϱ
Ϭ

ϱϬ

ϭϬϬ

ϭϱϬ

ϮϬϬ

FramesFig. 10 Profiles of iteration number required to satisfy theconvergence criterion of 0.01% for Fig. 9.

of viscoelastic fluid swells wider than the diameter ofan opening when the stream goes through the opening
because of viscoelastic forces which restore fluid’s original shape after the deformations at the opening. Although our velocity correcting method is based on approximated dynamics of viscoelastic fluid, we can plausibly generate Barus eﬀect.
5 Discussions and Limitations

Fig. 11 Spheres with diﬀerent viscoelasticity values successively thrown toward a solid dragon. Particles are color coded
based on their velocity correction coeﬃcients c (low coeﬃcient: cyan and high coeﬃcient: dark cyan).

Fig. 12 A viscoelastic ball thrown toward a solid cube.

4.3 Interactions
Fig. 11 illustrates several spheres with diﬀerent viscoelasticity values successively thrown toward a static
solid dragon. Similar to Figs. 4 (top) and 5 (right), ourmethod can easily handle merging of fluids and adhesion of fluids to a solid by updating pairwise interparticle connections.
Fig. 12 demonstrates two-way interactions of a viscoelastic material and a solid. The solid cube is gradually pulled toward a wall by the viscoelastic material
that connects to both of the solid cube and the wall.
This eﬀect cannot be simulated only with the previously proposed boundary-handling scheme for pressureand viscosity [1].

4.4 Barus Eﬀect
Fig. 13 demonstrates Barus eﬀect of a viscoelastic ball.
Barus eﬀect (also known as die swell, exclude swell,and Merrington eﬀect) is a phenomenon that a stream

Our method combines PBD with SPH to take advantage of geometric methods, such as numerical stability
and capability of handling large deformations. Additionally, unlike SPH-based methods, our method has
another benefit that time steps are not generally restricted by viscoelasticity values, namely velocity correction coeﬃcient c in our method; and thus our method
can use larger time steps and be faster than SPH-based
methods. In contrast to positive aspects, however, there
are a few negative issues on the adaptation of PBD.
First, our method exploits PBD to approximate viscoelastic behaviors instead of using physically-based models. Because of this approximation, our simulation results are less accurate than SPH-based methods which
use a physically-based model for viscoelastic eﬀects.
Second, we introduce several parameters which lack physical meanings into our method to cover various eﬀects
of viscoelastic fluids. These parameters need to be adjusted depending on each scenario through experiments
to generate desirable behaviors of materials. Third, since
results of our position-based velocity correction can beaﬀected by time steps similar to other geometric methods, our method would produce diﬀerent behavior ofviscoelastic fluids depending on time steps.
When materials are driven to extreme situations(e.g., a ball forcibly and significantly deformed by aheavy object), particle penetrations can be observed because of interparticle connections which exert larger velocity changes for farther connected particles to enablelarge deformations. Although this can be addressed byadjusting the influence of velocity corrections, such an
adjustment could make viscoelastic materials fail to restore their original shapes from highly deformed shapes.
Most viscoelastic materials undergo phase and property changes, e.g., turning into a stiﬀ solid, as timeVolume Preserving Viscoelastic Fluids with Large Deformations Using Position-based Velocity Corrections

passes. To simulate stiﬀ materials, strengthening theinfluence of velocity corrections would cause instability
in simulations. Thus, investigating how to simulate softand hard materials with phase and property changes ina unified framework is our important future work.
References
1. Akinci, N., Ihmsen, M., Akinci, G., Solenthaler, B.,
Teschner, M.: Versatile rigid-fluid coupling for incompressible SPH. ACM Transactions on Graphics 31(4),62:1–62:8 (2012)
2. Becker, M., Ihmsen, M., Teschner, M.: Corotated SPH
for deformable solids. In: Proceedings of the Fifth Eurographics Conference on Natural Phenomena, pp. 27–34(2009)
3. Becker, M., Teschner, M.: Weakly compressible SPH for
free surface flows. In: Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 209–217 (2007)
4. Bender, J., Müller, M., Otaduy, M.A., Teschner, M.:
Position-based methods for the simulation of solid objects in computer graphics. In: EUROGRAPHICS 2013
State of the Art Reports, pp. 1–22 (2013)
5. Bodin, K., Lacoursiere, C., Servin, M.: Constraint fluids. IEEE Transactions on Visualization and Computer
Graphics 18(3), 516–526 (2012)
6. Chang, Y., Bao, K., Liu, Y., Zhu, J., Wu, E.: A particlebased method for viscoelastic fluids animation. In: Proceedings of the 16th ACM Symposium on Virtual Reality
Software and Technology, pp. 111–117 (2009)
7. Clavet, S., Beaudoin, P., Poulin, P.: Particle-based viscoelastic fluid simulation. In: Proceedings of the 2005
ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 219–228 (2005)
8. Cummins, S.J., Rudman, M.: An SPH projection
method. Journal of Computational Physics 152(2), 584–
607 (1999)
9. Gerszewski, D., Bhattacharya, H., Bargteil, A.W.: A
point-based method for animating elastoplastic solids. In:
Proceedings of the 2009 ACM SIGGRAPH/Eurographics
Symposium on Computer Animation, pp. 133–138 (2009)
10. He, X., Liu, N., Li, S., Wang, H., Wang, G.: Local poisson
SPH for viscous incompressible fluids. Computer Graphics Forum 31(6), 1948–1958 (2012)
11. Ihmsen, M., Akinci, N., Becker, M., Teschner, M.: A parallel SPH implementation on multi-core CPUs. Computer
Graphics Forum 30(1), 99–112 (2011)
12. Ihmsen, M., Cornelis, J., Solenthaler, B., Horvath, C.,
Teschner, M.: Implicit incompressible SPH. IEEE Transactions on Visualization and Computer Graphics 20(3),426–435 (2014)
13. Ihmsen, M., Orthmann, J., Solenthaler, B., Kolb, A.,
Teschner, M.: SPH fluids in computer graphics. In: EUROGRAPHICS 2014 State of the Art Reports, pp. 21–42(2014)
14. Jones, B., Ward, S., Jallepalli, A., Perenia, J., Bargteil,
A.W.: Deformation embedding for point-based elastoplastic simulation. ACM Transactions on Graphics 33(2),21:1–21:9 (2014)
15. Macklin, M., Müller, M.: Position based fluids. ACM
Transactions on Graphics 32(4), 104:1–104:5 (2013)
16. Macklin, M., Müller, M., Chentanez, N., Kim, T.Y.: Unified particle physics for real-time applications. ACM
Transactions on Graphics 33(4), 104 (2014)917. Mao, H., Yang, Y.H.: Particle-based non-Newtonian fluid
animation with heating eﬀects. Tech. rep. (2006)
18. Miller, G., Pearce, A.: Globular dynamics: A connectedparticle system for animating viscous fluids. Computers& Graphics 13(3), 305 – 309 (1989)
19. Müller, M., Charypar, D., Gross, M.: Particle-based fluid
simulation for interactive applications. In: Proceedingsof the 2003 ACM SIGGRAPH/Eurographics Symposiumon Computer Animation, pp. 154–159 (2003)
20. Müller, M., Heidelberger, B., Hennix, M., Ratcliﬀ, J.: Position based dynamics. Journal of Visual Communication
Image Representation 18(2), 109–118 (2007)
21. Müller, M., Heidelberger, B., Teschner, M., Gross, M.:
Meshless deformations based on shape matching. ACM
Transactions on Graphics 24(3), 471–478 (2005)
22. Müller, M., Keiser, R., Nealen, A., Pauly, M., Gross, M.,
Alexa, M.: Point based animation of elastic, plastic andmelting objects. In: Proceedings of the 2004 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 141–151 (2004)
23. Solenthaler, B., Pajarola, R.: Predictive-corrective incompressible SPH. ACM Transactions on Graphics
28(3), 40:1–40:6 (2009)24. Solenthaler, B., Schläfli, J., Pajarola, R.: A unified particle model for fluid-solid interactions. Computer Animation and Virtual Worlds 18(1), 69–82 (2007)
25. Steele, K., Cline, D., Egbert, P.K., Dinerstein, J.: Modeling and rendering viscous liquids. Computer Animation
and Virtual Worlds 15(3-4), 183–192 (2004)
26. Takahashi, T., Fujishiro, I., Nishita, T.: A velocity correcting method for volume preserving viscoelastic fluids.
In: Proceedings of the Computer Graphics International
2014 (2014)
27. Takahashi, T., Nishita, T., Fujishiro, I.: Fast simulationof viscous fluids with elasticity and thermal conductivity
using position-based dynamics. Computers & Graphics
43(0), 21 – 30 (2014)
28. Takamatsu, K., Kanai, T.: A fast and practical method
for animating particle-based viscoelastic fluids. The International Journal of Virtual Reality 10, 29–35 (2011)
29. Tamura, N., Nakaguchi, T., Tsumura, N., Miyake, Y.:
Spring-bead animation of viscoelastic materials. IEEE
Computer Graphics and Applications 27(6), 87–93(2007)
30. Terzopoulos, D., Platt, J., Fleischer, K.: Heating andmelting deformable models. Journal of Visualization andComputer Animation 2, 68–73 (1991)
31. Zhou, Y., Lun, Z., Kalogerakis, E., Wang, R.: Implicit
integration for particle-based simulation of elasto-plastic
solids. Computer Graphics Forum 32(7), 215–223 (2013)