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

A Velocity Correcting Method for Volume Preserving
Viscoelastic Fluids
Tetsuya Takahashi · Issei Fujishiro · Tomoyuki Nishita

Abstract We propose a new particle-based method for
simulating viscoelastic fluids which preserve their volumes. Our method achieves the volume preservation
by enforcing the incompressibility of fluid, while correcting particle velocities to approximate the dynamics
of viscoelastic fluids without disturbing computations
for the incompressible flow. We oﬀer three schemes forcorrecting particle velocities. The first scheme employs
Shape Matching proposed by Müller et al. to derive appropriate transformations of particle sets. The second
computes attraction forces on the basis of Hooke’s law
to restrict particle motions. The third utilizes PositionBased Dynamics to restore the original relations of particle positions. The first scheme enables smooth transfers of deformation waves, the second is intuitive andsimple, and the third is easy to tune parameters. We
demonstrate that our method can preserve fluid volumes while generating plausible viscoelastic motions.
Keywords Fluid simulation · Viscoelasticity · Volume
preservation · Velocity correction

1 Introduction
Viscoelastic fluids are common in our everyday lives,
as we often see examples of these materials: egg white,
Tetsuya Takahashi
UEI Research / Keio University
E-mail: tetsuya.takahashi@uei.co.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

gels, mucus, gelatin, and toothpaste. Due to the ubiquity of viscoelastic materials, simulating their sticky
and bouncy behaviors has been required for featurefilms and video games.
Preserving fluid volumes is one of the most important aspects to simulate plausible fluid motions. For example, a lump of viscoelastic fluid on a flat solid surface
when being compressed does not horizontally spread
unless the volume of the lump is suﬃciently preserved.
In the computer graphics community, there are several
known simulation methods for viscoelastic fluids which
preserve fluid volumes using Eulerian grids [14] or Lagrangian meshes [2, 42, 41, 9], yet we herein focus on the
Lagrangian particle-based approach due to its conceptual simplicity.
In the literature of particle-based methods, combinations of Smoothed Particle Hydrodynamics (SPH)[25] and a geometrically-motivated method have been
proposed to simulate viscoelastic fluids eﬃciently androbustly [10, 37]. However, these previous combinations
spoil the fast convergence rate of state-of-the-art fluid
solvers, e.g., Predictive-Corrective Incompressible SPH[32] (PCISPH) and Implicit Incompressible SPH [17](IISPH), which can preserve fluid volumes well by enforcing the incompressibility of fluid. Takamatsu andKanai [37] combined SPH with a variant of Shape Matching [27] (SM), and interpolated two velocities derived
from both methods to fast and robustly simulate viscoelastic fluids. However, since the variant of SM does
not preserve the object volumes, the combined methodin [37] suﬀers from the loss of fluid volumes. A predictionrelaxation scheme proposed by Clavet et al. [10] directly
corrects particle positions based on various factors suchas collision responses and spring-based forces, yet thisscheme changes sets of neighboring particles and interparticle distances. Thus, their prediction-relaxation

2

Tetsuya Takahashi et al.

Fig. 1 A sphere of viscoelastic fluid compressed by a moving solid along a fixed path. Horizontally spreading motions areachieved due to the small volume loss of the material less than 1.0% by making use of position-based velocity corrections.

scheme disturbs the computations for the incompressible flow.
In this paper, we therefore propose a new particlebased method for simulating viscoelastic fluids which
preserve their volumes. Our method utilizes velocity
corrections to separate formulations for viscoelastic effects from computations for incompressible fluids, andour main contribution lies in separating the formulations not to slow down the computations for fluid incompressibility. To show our idea is applicable to diﬀerent geometrically-motivated methods and to compare
the existing methods with ours under fair conditions,
we oﬀer three schemes for correcting velocities of particles without estimating viscoelastic force fields for fastand stable simulations. The first scheme adopts SM[27] to derive suitable transformations of particle sets.
The second is based on Hooke’s law to compute attraction forces. The last utilizes the idea of Position-Based
Dynamics [26] (PBD) to compute velocity correction
vectors. We employ IISPH [17] as our underlying fluid
solver owing to the fast convergence of density fluctuations for the incompressible flow, and combine IISPH[17] with one of the three velocity correction schemes
to achieve volume preserving viscoelastic fluids. Fig. 1
illustrates horizontally spreading motions of a viscoelastic ball compressed by a moving solid, which cannot be
fully captured with compressible viscoelastic fluids.

2 Related Work
Viscoelastic fluids Viscoelastic fluids have mainly been
simulated using a spring-based method or methods basedon SPH. Miller and Pearce [22] and Terzopoulos et al.
[39] proposed a spring-based model that computes repulsion and attraction forces among particles to simulate viscoelastic materials. Their method was also adopted by Steele et al. [34] and Tamura et al. [38]. Clavet etal. [10] combined this spring-based method with SPH

to simulate materials with elasticity, plasticity, and viscosity, adopting a prediction-relaxation scheme, which
partly shares the features of PBD. Takamatsu and Kanai[37] fast and robustly simulated viscoelastic materials
by combining SPH with a variant of SM. Müller et al.
[28] proposed an elastic term which uses Moving Least
Square (MLS) to simulate elastoplastic objects. Solenthaler et al. [33] adopted the formulation of this elastic
term and computed it using SPH instead of MLS toallow for robustly simulating fluid with various properties under the condition of collinear or coplanar particle
distributions. The method of Solenthaler et al. [33] wasextended to handle rotational motions of elastic materials [3]. Mao et al. [21] introduced an elastic force term,
called nonlinear corotational Maxwell model, into theNavier-Stokes equations to simulate viscoelastic fluids.
This method was also adopted by Chang et al. [8]. Paiva
et al. [29] simulated viscoplastic objects using a generalized Newtonian model.
Gerszewski et al. [13] 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. [43]in an implicit manner to eﬃciently and robustly perform simulations adopting larger time steps.
Takahashi and Fujishiro [35] simulated viscous fluidswith elasticity using PBD with a low computational
cost.
Volume preservation Volume preservation is necessary to plausibly generate motions of objects, andmany methods for preserving volumes have been proposed for deformable objects by computing deformationenergy [40, 30] or imposing constraints upon computational elements [16, 19, 36, 12].
For particle-based fluids, volume preservation methods have focused on minimizing density fluctuations
from the fluid rest density by enforcing the incompressibility of fluid [18], because of frequent topology changes
unlike methods for deformable objects. After Müller etA Velocity Correcting Method for Volume Preserving Viscoelastic Fluids

al. [25] presented the basic SPH method, which suﬀers
from unacceptable density fluctuations and volume loss
due to the dependence on an equation of state (EOS),
Becker and Teschner [4] proposed Weakly Compressible SPH and alleviated the fluctuations by replacing
the EOS in [25] with a stiﬀer EOS called Tait equation[24], necessitating excessively smaller time steps. In order to adopt larger ones, predictive-correction schemes
were proposed, e.g., PCISPH [32] and local Poisson
SPH [15]. Recently, Bodin et al. [6] proposed a system of velocity constraints to improve the accuracy of
PCISPH achieving uniform density fields over particles.
To further accelerate fluid simulations using PCISPH,
Macklin and Müller [20] adapted PBD for fluid simulation to allow for the use of larger time steps than theones adopted in PCISPH. Similarly, Ihmsen et al. [17]proposed IISPH, which computes particle pressure values in an implicit manner with minimal assumptions.
Furthermore, Cornelis et al. [11] combined IISPH with
Fluid-Implicit-Particle to further accelerate simulation
performances.

3 Our Algorithm

3

where vi is the velocity of particle i and F̃i is nonpressure force, using the continuity equation by∑ρ̃i = ρi + ∆t
mj ṽij ∇Wij ,jwhere ∆t is a time step, and ṽij = ṽi − ṽj . The next
step is to find pressure values which result in appropriate pressure forces to make all particles satisfy the rest
density ρ0 , and we obtain the following equations:ρ0 − ρ̃i = pi aii + bij ,∑aii =mj (dii − dji ) · ∇Wij ,
j

bij =

∑

∑
∑
mj (
dij pj − djj pj −djk pk ) · ∇Wij ,jdii = −∆t

j

2

∑ mj
j

dij = −∆t2

ρ2i

k̸=i

∇Wij ,

mj∇Wij .
ρ2j

Then, we compute pressure values using weighted Jacobi method with an iteration index n and a relaxation
factor w:
ρ0 − ρ̃i − bnij
pn+1
= (1 − w)pni + w.
iaii

Our key idea to simulate volume preserving viscoelastic
fluids is to separately deal with the volume preservationand viscoelasticity. While preserving the fluid volumes
by enforcing the incompressibility of fluid using IISPH[17], we utilize velocity corrections for viscoelastic effects without negatively influencing the convergence ofdensity fluctuations in IISPH [17].
In Section 3.1, we first make a brief explanationof IISPH [17] to elucidate the mechanism of the fluid
solver and its important features, and then Section 3.2
shows the procedure of our method. The three velocity
correction schemes are detailed in Section 4.

Note that the computations of the weighted Jacobi method include particle positions, and thus the positions
must be kept during the iterations after intermediate
velocity ṽ is determined. As explained in [32] and [17],ignoring position changes of particles is the cause oferroneous sets of neighboring particles and their interparticle distances, and leads to inaccurate estimation ofphysical values, which can cause the fluid solver to failor slow down the convergence of the iterations.
After we obtain pressure pi , pressure force Fpi iscomputed by∑pjpi
Fpi = −mi(1)
mj ( 2 + 2 )∇Wij ,
ρ
ρ
i
j
j

3.1 IISPH Method

and particle velocity vi and position xi are integratedusing the Euler-Cromer scheme.
Since we get pressure p from intermediate velocity ṽ,IISPH is interpreted as a function that p = IISP H(ṽ),
and this corresponds to the pressure projection used in
the Eulerian fluid simulation [7].

IISPH aims to achieve deviations of the particle density
ρ less than a certain percentage, e.g. 1%, to the fluid
rest density ρ0 by iteratively correcting the pressure ofeach particle in an implicit manner. For the derivationof formulations and implementation details, please referto the original paper of IISPH [17].
∑First, we compute the particle density ρi = j mj Wij
for particle i with its mass mi , neighboring particle j,
and finite support kernel Wij = W (xij , h), where xij =xi − xj , xi is the position of particle i, and h denotes
a kernel radius. Then, intermediate density ρ̃i is computed with intermediate velocity ṽi = vi + ∆tF̃i /mi ,3.2 Procedure of Our Method
We insert the velocity correction step (Section 4) before the iterations in IISPH, and the following is our
simulation algorithm.
1. Find neighboring particles;4

2. Correct particle velocity (see Section 4);(a) Compute predicted velocity vadv ;(b) Get corrected velocity v∗ using one of the threeschemes;(c) Obtain intermediate velocity ṽ using XSPH [23];3. Compute pressure value p = IISP H(ṽ);
4. Compute pressure force Fp using Eq. (1);
5. Update velocity v = ṽ + ∆tFp /m;
6. Update position x = x + ∆tv.

4 Velocity Correction
We correct particle velocities to describe viscoelastic
eﬀects without changing particle positions. In order toobtain intermediate velocity ṽi for particle i as an input
for IISPH, we first compute predicted velocity viadv withall forces Fadv
except for viscoelastic and pressure onesi
without directly computing ṽi as in IISPH:
Fadv
i
.
mi
Then, we correct particle velocities with velocity correction vector ∆vi :

viadv = vi + ∆tvi∗ = viadv + ∆vi .
Lastly, we apply XSPH [23] to reduce particle oscilla∗= vi∗ − vj∗ and a velocity mixing paramtions with vij
eter ϵ (0 ≤ ϵ ≤ 1):∑ mj∗
ṽi = vi∗ + ϵ
vji
Wij .
ρ
j
j
Note that XSPH [23] modifies only particle velocities
without aﬀecting the convergence of the fluid solver,
and therefore XSPH [23] is suitable for the purpose of
our velocity correction.
In the followings, we explain the three schemes forcomputing velocity correction vector ∆v.

4.1 Shape Matching Scheme
We formulate shape matching scheme as with the methodof [37] adapting Lattice Shape Matching [31], which was
extended from the original SM [27] to increase the degree of freedom of object deformations.
We consider particle i with its reference position
x0i , which has a set of Sism particles for SM, and thereference position s0i of the set is defined to be the center of mass of the particles. In SM, particle i is pulled
toward its goal position gi to restore the original configuration of the particles, and individual goal positions
are computed to match the original configuration of theparticles defined by x0i with current particle distributions denoted as xi after the particles are transferred.

Tetsuya Takahashi et al.

Specifically, we compute a rotational matrix Ri whichtransforms the particles to their own reference positions
with si , the current center of mass of the particles. The
goal position gi is computed by
S sm

i
(
)
1 ∑
Rj (x0i − x0j ) + sj
gi = smSij= Ri (x0i − s0i ) + si ,and then we obtain velocity correction vector ∆vism
with a coeﬃcient for the shape matching velocity correction kism (0 ≤ lism ≤ kism ≤ 1) and lism , the lower
limit of kism :gi − x i,∆t
If lism < kism , and αism < (||gi − xi ||)/h, where αism (0 ≤αism ) is a yield criterion, which invalidates the influenceof small oscillations and deformations of materials, we
weaken the eﬀect of kism to realize plastic and viscouseﬀects by∆vism = kism

sm
kism ← max(kism − ∆tkism dsm
i , li ),where dsm
(0 ≤ dsmii ) is a parameter that controls theweakening speed of kism .
In order to enable the merge of viscoelastic materials, we change a set of particles for SM. Let csm
ij denotethe distance of two particles which are not associatedwith each other as a particle for SM, and we updatesmsmsm
the set if csm+ βjsm )/2,
ij /h < βij , where βij = (βi(0 < βism ), is a threshold to update a set of particles
for SM. We also update the set for the split of the materials if γism < ||gi − xi ||/h, where γism (0 < γism ) is aparameter which limits the extent of deformations.

4.2 Spring-based Scheme
In the spring-based scheme, we correct particle velocities with springs between two particles, and the springs
are created when an object is generated, Hereafter wecall particles connected with their springs as connectedparticles.
Since IISPH that we use can resolve fluid compression, we address only fluid expansions which are causedby separating particles from the other. In order to compute velocity correction vector vispr only with connectedparticles which are separating from the other, we definea distance function with their positions and interparticle distance rij as
D(xi , xj ) = max(||xij || − rij , 0),(2)where rij is initialized with the distance between thetwo particles when their spring is created. Then, we

A Velocity Correcting Method for Volume Preserving Viscoelastic Fluids

compute velocity correction vector vispr based on Hooke’s
law with kispr (0 ≤ lispr ≤ kispr ), a coeﬃcient for thespring-based velocity correction and lispr , the lower limitof kispr by∆t∆vispr = −
mi

Sispr

∑ kispr + kjspr
2j

D(xi , xj )xij,||xij ||

where Sispr denotes the number of connected particles
to particle i. Similar to the shape matching scheme, we
weaken the eﬀect of kispr if lispr < kispr :5

Table 1 Simulation conditions and performances.
Fig.
1
2 (a)
2 (b)
2 (c)
3 (a)
3 (b)
5 (a)
5 (b)
7
8

# of particles
fluid/solid
65.8k/45.0k
17.1k/23.3k
17.1k/23.3k
17.1k/23.3k
8.2k/13.4k
8.2k/13.4k
8.2k/13.4k
8.2k/13.4k
33.4k/40.7k
45.2k/36.4k

Total
time
10m:16s
5m:53s
4m:51s
4m:55s
6m:28s
6m:46s
6m:58s
6m:31s
12m:15s
58m:05s

Total
simulation steps
60
150
150
150
200
200
150
150
200
170

S spr

i
∆t ∑
spr
kispr ← max(kispr − spr
dspr
ij σij , li ),Sij{D(xi ,xj )spr
1 if αij
<
,
h
σij =
0 otherwise,

spr
spr
where dspr
+ dsprij = (dij )/2, (0 ≤ di ), governs thespr
spr
weakening speed of ki , and αij = (αispr + αjspr )/2,(0 ≤ αispr ), is a yield criterion.
In the spring-based scheme, we also enable the merge
and split of viscoelastic materials. Let cspr
ij denote thedistance of two particles which are not interconnected,
and we generate a new spring between the two parspr
spr
spr
ticles if cspr
+ βjspr )/2,
ij /h < βij , where βij = (βispr(0 < βi ), is a threshold to generate springs. By contrast, we remove the spring between two particles if
spr
spr
γij
< ||xij − rij ||/h, where γij= (γispr + γjspr )/2,spr
(0 < γi ), is a parameter which restricts the extent of
deformations.

4.3 Position-based Scheme
The position-based scheme is similar to the spring-based
one since particle velocities are corrected based on pairwise connections between particles, and the connections
are added and removed as in the spring-based scheme.
The diﬀerences between the two schemes lie in how tocorrect particle velocities and select relevant parameters. This position-based scheme relies on positional
diﬀerences of particles to compute velocity correction
vector ∆vipb by adapting the idea of PBD [26,5] instead of computing attraction forces:pb

1∆vipb = −∆tSi∑
kipb + kjpb
j

2

mj
xij
D(xi , xj ),mi + mj||xij ||

where Sipb is the number of connected particles to particle i, kipb (0 ≤ lipb ≤ kipb ≤ 1) is a coeﬃcient for the
position-based velocity correction, and lipb , the lower
limit of kipb . Since the range of kipb is bounded by twosides, tuning parameters can be easier in the positionbased scheme than the spring-based one. We correct

particle velocities and modify k pb as in the spring-based
scheme.

5 Results
We implemented our method in C++ and parallelized itwith OpenMP 2.0. All the simulations were performed
on a PC with a 4-core Intel Core i7 3.50 GHz CPU
and RAM 16.0 GB. We gained accelerated simulation
performances by a factor from three to five with six
threads in total. In our simulations, fluid-solid couplingis achieved by the method of Akinci et al. [1]. For fast
simulations, we used only 2 iterations for computing
pressure values in IISPH except Figs. 1, 3, and 5. In all
the simulations, parameters were empirically adjusted.
Our velocity correction schemes occupy roughly 16.4%of the total computational time with 2 pressure iteration in IISPH. We used oﬀ-line renderer POV-Ray 3.7
and the rendering needed about 20 seconds per frameon average for Fig. 1. The accompanying video includes
several results produced with our method or existing
methods. We tabulate the simulation conditions andperformances in Table 1.

5.1 Velocity Correction Schemes
We compare the motions of viscoelastic fluids producedwith the three diﬀerent schemes for velocity corrections
in order to clarify how the selection of the schemes affects the resulting behaviors. Fig. 2 displays the motionof a viscoelastic ball dropped onto the ground, where
we used the same number of neighboring particles in all
the schemes, and we can plausibly generate characteristic behaviors of viscoelastic fluids. Our method using
the shape matching scheme (Fig. 2 (a)) enables smoothtransfers of deformation waves, which leads to preserving energy for a longer time as compared to the others.
Since the three schemes are not physically-based, we6Tetsuya Takahashi et al.

(a) Shape matching scheme(a) Method of Takamatsu and Kanai [37](b) Spring-based scheme(b) Our method using the shape matching scheme(c) Position-based scheme
Fig. 2 Motion comparisons of the three diﬀerent velocity
correction schemes.

need to carefully adjust simulation parameters to generate similar results with diﬀerent schemes.
As shown in Table 1, the simulation costs of themethods using the spring-based scheme (Fig. 2 (b)) andthe position-based one (Fig. 2 (c)) are almost the same.
However, the method using the shape matching scheme
can be a bit costly due to the singular value decomposition in SM, and besides needs more particles in theset for SM to stabilize the simulation than the others.

5.2 Volume Preservation
We compare the existing methods and ours throughinteractions between a viscoelastic ball with a moving
solid along a fixed path.
Shape matching scheme First, we elucidate thediﬀerence between the method of Takamatsu and Kanai[37], who used the combination of SPH and SM, andours using the shape matching scheme for velocity corrections. Although Takamatsu and Kanai [37] employed
the basic SPH as their underlying fluid solver, which
suﬀers from significant volume loss, we replace the SPH
with IISPH for a fair comparison. Namely, we here regard the method in [37] as the combination of IISPH

Fig. 3 A volume comparison with a viscoelastic sphere compressed by a moving solid along a fixed path.

and SM. In this scene, we used 5 pressure iterations
in IISPH, and we estimated an occupied volume
V
∑
V
by
a
summationofparticles’
volume
as
V
=
i i =
∑
m
/
ρ
ˆ
,
where
V
is
the
volumeofparticle
i,
and
ρˆi
i
i
i
i
the estimated density after only the particle i is moved
using Laplacian smoothing to reduce the eﬀect of theunderestimated particle density on fluid surfaces.
Fig. 3 compares the two methods, and Fig. 4 shows
the profile of the comparison in terms of the volume
V . With the method of Takamatsu and Kanai [37], thevolume of the ball fluctuates due to the first collision
with the ground, and decreases and increases according
to the movement of the solid. In contrast, our method
can preserve the volume of the ball, which is fairly close
to the reference value, and achieve less fluctuations in
the volume. Additionally, the ball exhibits horizontally
spreading motions when being compressed, which cannot be captured fully with the method in [37].
Spring-based scheme Second, we compare themethod of Clavet et al. [10], who combined SPH with aspring-based method adopting a prediction-relaxation
scheme, with ours using the spring-based scheme. As
with the above previous method, since Clavet et al. [10]used the basic SPH, we again use IISPH as the underlying fluid solver of their method. In this scenario, we
set the volume loss of 0.01% as a convergence criteriain IISPH to clarify the diﬀerences of the two methods.

A Velocity Correcting Method for Volume Preserving Viscoelastic Fluids

7

200

0.055

180
160

0.045

140
0.035

Takamatsu and Kanai
0.025

Ours
Reference
20

40

60

80

100

120

140

160

120
100
80

0.015
0

Iterations

Volume V (m3)

220
0.065

180

60

Simulation steps

Fig. 4 Changes in the volume of the viscoelastic sphere in
Fig. 3 according to simulation steps.

40
Clavet et al.
20

Ours

0
0

20

40

60

80

100

Simulation steps

Fig. 6 A profile of the number of needed pressure iterations
to achieve the volume loss of the material less than 0.01% in
Fig. 5 according to simulation steps.

(a) Method of Clavet et al. [10]Fig. 7 A ball with diﬀerent viscoelasticity values dropped
onto a funnel (low viscoelascitiy: light purple, and high one:dark purple).

5.3 Diﬀerent Viscoelasticity Values(b) Our method using the spring-based scheme
Fig. 5 A comparison with a viscoelastic sphere compressedby a moving solid along a fixed path.

Fig. 5 compares the two methods, and Fig. 6 gives
the profile of the number of the needed iterations in
Fig. 5. The method of Clavet et al. [10] changes particle positions in the prediction-relaxation scheme, andthese changes make the fluid solver fail or slow down
convergence of the iterations. As illustrated in Fig. 6,
the method of Clavet et al. [10] requires more pressure
iterations to enforce the incompressibility of the material than our method when the material is compressedby the solid, because spring forces tend to be strong andlead to significant position displacements. Additionally,
though IISPH can produce satisfactory pressure values
to achieve the incompressibility of fluid, resulting positions through the time integration are not optimal, andhence can cause significant oscillations in [10].

Fig. 7 demonstrates the eﬀect of diﬀerent viscoelasticity values. In this scene, a ball is dropped onto a funnel. While the balls with lower viscoelasticity valuessmoothly fall through the funnel, the highly viscoelastic ball needs longer time to pass through the funnel.
Additionally, in Figs. 7 (b) and (c), we can observe theBarus eﬀect (also known as die swell, exclude swell, andMerrington eﬀect) that a stream of viscoelastic fluids
swells wider than the diameter of an opening when thestream goes through the opening. Our velocity correction method can plausibly generate such a special eﬀectof viscoelastic fluids.
Fig. 8 illustrates the merge of the two balls with different viscoelasticity values. Due to the update of particle connections, we can easily and naturally generateadhesive behavior of the two viscoelastic balls.

6 Conclusions
We have presented a new particle-based method for
simulating volume preserving viscoelastic fluids by correcting particle velocities based on particle positions.
The volume preservation is achieved by separating for-

8

Tetsuya Takahashi et al.

Fig. 8 Two balls with diﬀerent viscoelasticity values colliding with each other (low viscoelasticity: light purple and highone: dark purple).

mulations for viscoelastic eﬀects from computations forincompressible fluids. To correct particle velocities, wehave provided three diﬀerent schemes: shape matching,spring-based, and position-based velocity corrections.
We have empirically proven that the three schemes canplausibly simulate viscoelastic fluids, and that our method can preserve fluid volumes generating characteristic
viscoelastic motions such as Barus eﬀect.
Since our method utilizes the geometric relationsof particles to correct particle velocities with several
simulation parameters instead of estimating viscoelastic force fields, we cannot adopt real physical properties,
and needs to empirically adjust the parameters throughexperiments in order to generate desirable viscoelastic
eﬀects. For future work, we plan to address the problem by the control of fluid behavior, and investigating
further how to find appropriate simulation parameters.

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),
1–8 (2012)
2. Bargteil, A.W., Wojtan, C., Hodgins, J.K., Turk, G.: A finite element method for animating large viscoplastic flow.
ACM Transactions on Graphics 26(3) (2007)
3. 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)
4. 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)
5. Bender, J., Müller, M., Otaduy, M.A., Teschner, M.:
Position-based methods for the simulation of solid ob-

jects in computer graphics. In: EUROGRAPHICS 2013
State of the Art Reports, pp. 1–22 (2013)
6. Bodin, K., Lacoursiere, C., Servin, M.: Constraint fluids. IEEE Transactions on Visualization and Computer
Graphics 18(3), 516–526 (2012)
7. Bridson, R.: Fluid simulation for computer graphics. A
K Peters / CRC Press (2008)
8. 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)
9. Clausen, P., Wicke, M., Shewchuk, J.R., O’Brien, J.F.:
Simulating liquids and solid-liquid interactions with lagrangian meshes. ACM Transactions on Graphics 32(2),1–15 (2013)
10. 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)
11. Cornelis, J., Ihmsen, M., Peer, A., Teschner, M.: IISPHFLIP for incompressible fluids. Computer Graphics Forum 33(2), 255–262 (2014)12. Diziol, R., Bender, J., Bayer, D.: Robust real-time deformation of incompressible surface meshes. In: Proceedingsof the 2011 ACM SIGGRAPH/Eurographics Symposiumon Computer Animation, pp. 237–246 (2011)
13. 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)
14. Goktekin, T.G., Bargteil, A.W., O’Brien, J.F.: A method
for animating viscoelastic fluids. ACM Transactions on
Graphics 23(3), 463–468 (2004)
15. 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)16. Hong, M., Jung, S., Choi, M.H., Welch, S.W.: Fast volume preservation for a mass-spring system. IEEE Computer Graphics and Applications 26(5), 83–91 (2006)
17. 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)
18. 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)
19. Irving, G., Schroeder, C., Fedkiw, R.: Volume conserving
finite element simulations of deformable models. ACM
Transactions on Graphics 26(3) (2007)
20. Macklin, M., Müller, M.: Position based fluids. ACM
Transactions on Graphics 32(4), 1–5 (2013)
21. Mao, H., Yang, Y.H.: Particle-based non-newtonian fluid
animation with heating eﬀects. Tech. rep., University ofAlberta (2006)
22. Miller, G., Pearce, A.: Globular dynamics: A connectedparticle system for animating viscous fluids. Computers
and Graphics 13(3), 305–309 (1989)
23. Monaghan, J.J.: On the problem of penetration in particle methods. Journal of Computational Physics 82(1),
1–15 (1989)
24. Monaghan, J.J.: Simulating free surface flows with
SPH. Journal of Computational Physics 110(2), 399–406(1994)
25. 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)

A Velocity Correcting Method for Volume Preserving Viscoelastic Fluids
26. Müller, M., Heidelberger, B., Hennix, M., Ratcliﬀ, J.: Position based dynamics. Journal of Visual Communication
and Image Representation 18(2), 109–118 (2007)
27. Müller, M., Heidelberger, B., Teschner, M., Gross, M.:
Meshless deformations based on shape matching. ACM
Transactions on Graphics 24(3), 471–478 (2005)
28. 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)
29. Paiva, A., Petronetto, F., Lewiner, T., Tavares,
G.: Particle-based viscoplastic fluid/solid simulation.
Computer-Aided Design 41(4), 306–314 (2009)30. Patterson, T., Mitchell, N., Sifakis, E.: Simulation ofcomplex nonlinear elastic bodies using lattice deformers. ACM Transactions on Graphics 31(6), 197:1–197:10(2012)
31. Rivers, A.R., James, D.L.: FastLSM: Fast lattice shape
matching for robust real-time deformation. ACM Transactions on Graphics 26(3) (2007)
32. Solenthaler, B., Pajarola, R.: Predictive-corrective incompressible SPH. ACM Transactions on Graphics
28(3), 1–6 (2009)33. 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)
34. 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)
35. Takahashi, T., Fujishiro, I.: Accelerated viscous fluid simulation using position-based constraints. In: Proceedingsof CAD/Graphics 2013, pp. 260–267 (2013)
36. Takamatsu, K., Kanai, T.: Volume-preserving LSM deformations. In: ACM SIGGRAPH ASIA 2009 Sketches,
pp. 15:1–15:1 (2009)
37. 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)
38. Tamura, N., Nakaguchi, T., Tsumura, N., Miyake, Y.:
Spring-bead animation of viscoelastic materials. IEEE
Computer Graphics and Applications 27(6), 87–93(2007)
39. Terzopoulos, D., Platt, J., Fleischer, K.: Heating andmelting deformable models. Journal of Visualization andComputer Animation 2, 68–73 (1991)
40. Teschner, M., Heidelberger, B., Muller, M., Gross, M.: A
versatile and robust model for geometrically complex deformable solids. In: Proceedings of the Computer Graphics International, pp. 312–319 (2004)
41. Wicke, M., Ritchie, D., Klingner, B.M., Burke, S.,
Shewchuk, J.R., O’Brien, J.F.: Dynamic local remeshing for elastoplastic simulation. ACM Transactions on
Graphics 29, 1–11 (2010)
42. Wojtan, C., Turk, G.: Fast viscoelastic behavior withthin features. ACM Transactions on Graphics 27(3), 1–8(2008)
43. 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)9