7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Modeling and simulation of viscoelastic fluids using Smoothed Particle Hydrodynamics
Robert Rundqvist, Andreas Mark, Fredrik Edelvik and Johan S Carlson
Fraunhofer Chalmers Research Centre for Industrial Mathematics
Chalmers Science Park, SE412 88 G6teborg, Sweden
robert.rundqvist@fcc.chalmers.se
Keywords: Viscoelastic fluids, Thixotropy, Smoothed Particle Hydrodynamics, Free surface flow, Car body sealing
Abstract
Multiphase flow simulation using Smoothed Particle Hydrodynamics (SPH) has gained interest during recent years, mostly
due to the inherent flexibility of the method and the physically rather intuitive formulation of extra constitutive equations
needed when dealing with for instance nonNewtonian flows. In the work presented here, simulations based on an SPH model
implemented in the flow solver IBOFlow has been used for simulation of robotic application of sealing material on a car body.
Application of sealing materials is done in order to prevent water leakage into cavities of the body, and to reduce noise. In
offline programming of the robots in the automotive paintshop it is of great interest to predict shape and appearance of sealing
material without having to resort to trialanderror procedures.
The flow of sealing material in the air between applicator and target (car body) is relatively uncomplicated, as the material
mostly moves at constant velocity until impact on target. The flow of the material on the target is however more complex,
applied material flows at the target surface due to inertia, gravity and pressure and in order to predict the appearance of the
applied material, flow equations for a nonNewtonian fluid with a free surface needs to be solved. The sealing material is both
thixotropic and viscoelastic; the material is shear thinning but needs to be sheared for some time before the structure of the
material is broken down. Conversely, the regain of structure of the material, and thereby also the increase of viscosity when
shearing is stopped or reduced, is also connected to a delay time. In the model used, the local viscosity is considered obeying a
first order differential equation where the stationary limit is determined by a Bingham relation.
The simulation model was built by comparing simulations and experiments at three different stages of complexity. In the most
fundamental stage the material properties were determined. Using a rotational rheometer, yield stress, plastic viscosity and
thixotropy time constant were determined and implemented in the simulation model. To verify the numerical behaviour of the
rheology, simulated rheometer experiments were carried out and compared with the physical experiments. In the second stage,
simulation of application of sealing material with a stationary hollowcone nozzle was carried out. To verify the simulations,
the resulting thickness, width and shape of applied material as a function of time were compared to experiments. In the third
stage a moving applicator of the same type was considered, here thickness width and shape of applied material as a function of
applicator to target distances were compared between experiments and simulation. At all three stages the number of SPH
particles, i.e. grid points, in the simulations was varied in order to verify that the simulations were resolution independent
Results of the simulations show good agreement between experiments and simulations in all stages using no artificial tuning of
the models, that is, all parameters used in the models are based on physical considerations. Furthermore, simulation time on a
desktop computer indicate that computational power required for industrially relevant cases is not prohibitively large, for the
most complex cases in this work simulation time did not exceed six hours.
Introduction It is of great interest to get an a priori prediction of how
well the applied sealing material covers the intended
Sealing material is applied to car bodies in seams to protect regions; when commencing the production of new car
the car from moisture, water leakage or to dampen noice. models it is presently necessary to manually program the
The application, which is performed by robots, typically robots for new application paths in a timeconsuming
takes place after the bodyinwhite has been treated with trialanderror procedure. With a simulation tool able to
Phosfate and Electrocoat dips, but before the different paint make predictions of sealing material application and flow it
layers (Primer, Base coat, Top coat/Clear coat) are applied, would be possible both to program the robots beforehand
An example of a sealing seam on a bodyinwhite can be (offline programming) and also to be able to reduce
seen in Figure 1. material waste.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Nomenclature
Characteristic length scale (m)
Gravitational acceleration (m s2)
Kernel smoothing length (m)
SPH node mass (kg)
Total number of SPH nodes
Pressure (Pa)
SPH distance function (m)
Time step (s)
Time (s)
Characteristic time scale (s)
Velocity (m s')
Characteristic velocity (m s')
Kernel weight function (m3)
Position (m)
Figure 1: Sealing material as applied to a car body.
The sealing applicator investigated in this work is a so
called hollowcone nozzle, in which a conical curtain of
sealing material is sprayed onto the target. The material is
broken up into droplets with a size less than a millimetre in
the applicator or just downstream of it. The interaction
between material droplets and air is normally negligible; the
computational challenge starts as the droplets impact the
target, transferring its momentum to the string that is being
formed on the surface. The sealing material is
nonNewtonian, the viscosity decays with increasing shear,
and thixotropic meaning that there is a certain timelag
connected to the change of viscosity. This timelag is due to
the internal breakdown and buildup of the intermolecular
connections that gives the material its viscous structure. To
accurately model the thixotropic nature of the flow and to
simplify the tracking of the free surfaces a Lagrangean
method is preferred.
The simulations in this work are based on Smoothed Particle
Hydrodynamics, a fluidfollowing method that
automatically includes transport of history dependent fluid
properties and propagation of a fluid/surrounding
(liquid/air) interface. The principle behind the discretization
concept in SPH is that the function values in the PDEs are
expressed as weighted averages of the closest nodes using a
smoothing kernel. In this way a continuous function
approximating the solution is obtained on a discretized
network in which the computational nodes are advected
with the fluid velocity.
To provide the simulation with data on the physical
parameters of the used sealing material experiments were
preformed in a rotary rheometer. The material used is
Betaguard, a PVC sealer commonly used in car industry. To
quantify the inlet conditions from the applicator in use a test
series with a stationary applicator was performed and to
verify the simulations a number of test plates with straight
strings of material were produced.
The simulation models are currently implemented as a
demonstrator module in the FCC software IPS Virtual Paint,
in which the user can import geometry descriptions, define
robot paths and applicator operating conditions. Results in
terms of sealing material on the included geometries are
displayed as the simulation progresses.
Greek letters
0 Arbitrary field variable ()
S Rate of deformation (s')
pt Dynamic viscosity (Pa s)
p Density (kg m3)
r Shear stress (Pa)
Subscripts
* Intermediate value
Asymptotic value
i SPH node number
thix Thixotropic
yield Yield value
Experimental setup
The experiments were carried out at Volvo Car Corporation
in Torslanda, Giteborg. In the first set of experiments, the
rheology of the material used was determined. Using a
rotary rheometer with a parallel disc, relations between
applied torque and deformation was determined. As the
sealing material was expected to exhibit thixotropic
properties, the material was subjected to a series of constant
deformation of successively higher rate. When the
maximum rate of deformation, 450 reciprocal seconds, had
been reached, the same levels were measured once more but
going in the other direction thereby it becomes possible to
detect differences in timescale for internal buildup and
breakdown of material structure. The setting for shear rate
used in the experiments is shown in Figure 2.
Figure 2:
of time.
500 1000 1500 2000 2500 3000 3500
Time [s]
Rheometer deformation rate setting as a function
o0
0
A second set of experiments was set up in order to
characterize the sealing applicator studied, that is to
determine material flow rate and inlet velocities. The setup
is shown in action in Figure 3. Here sealing material was
applied to a flat plate using a fixed hollowcone nozzle. The
amount of sealing material and distance from applicator to
target was varied, and the buildup of material on the target
was recorded,. The velocity of the material leaving the
nozzle was estimated using a high speed camera and the
mass flow of material at a given applicator setting was
determined using a simple bucket and stopwatch approach.
In the default setting used, the flow rate of material was 28.5
g/s and the density of the sealing material at ambient
conditions was 1080 kg/m3.
In a third set of experiments the applicator was allowed to
move in a direction parallel to the target plate, as is shown
in Figure 4. The distance between applicator and target plate
was varied for different test runs, ranging from 20 to 60
millimeters, and the width of the applied strings of material
was measured to be used for comparison with the
simulations.
Figure 4: String application of sealing material onto a flat
plate in the image the applicator moves from right to left
with a constant velocity. Distance between nozzle and target
is 30 mm.
Numerical method
The SPH type of discretizations has gained in popularity the
last decade. A comprehensible review of the method is given
by Monaghan (2005). The SPH implementation in this paper
is mostly based on the work of Cummins (1999), but with
the extension to an implicit formulation of the momentum
equations in order to handle the high viscosity correctly and
to be able to take longer time steps. Field variables are in
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
principle evaluated in a point by averaging the values of the
surrounding nodes, weighted by the smoothing kernel. The
principle is illustrated in Figure 5 and the evaluation of a
variable 0 at a point x using a cubic spline kernel is shown
in equation (1).
0 o
I
Figure 5: SPH principles. A field variable is to be
evaluated at the position of the red SPH particle in the
figure. Using the smoothing kernel, the value will be
averaged using the nearest neighbours.
N M
O(X): =__ (x)w, (X)
71 P X
(2 q 4(1q)
,0 q 1
,l< q < 2
,q>2
q = xx I/h
Starting from NavierStoke's equations in the form given in
equation (2),
V.u=0
u UVu Vp+ V2 (2)
at p p
momentum and pressure equations are solved using a
predictorcorrector scheme. First, the neighbours of each
node is located through a kd tree search and local densities
pi(t) are calculated from the SPH node positions xi at the
time t. The intermediate velocities u*, of the SPH nodes are
then calculated by solving the momentum equations with
the pressure gradient set to zero:
ui _V 2i At = ui+gAt
IAi
Note that the discretization of this equation, specifically the
stress tensor which is expressed using the SPH formalism
outlined in equation (1), is handled implicitly in order to
handle the high viscosities in the system with a reasonable
timestep. From the solution of the linearized system we get
the intermediate velocity, which is used to update the node
positions and to calculate an intermediate density, p*i(t). The
next step is to find the pressure gradient of the next
timestep by solving equation (4).
V VP(t + At) =
p0*i oAt2
Again, the gradients in the equation are evaluated according
to SPH principles. Note that the unsubscripted p in
equation (4) is the real, incompressible, fluid density. In the
last step of the procedure the rest of the momentum equation
is solved, that is, the pressure gradient is compensated for:
Using a CrankNicholson scheme, positions at time t+At can
now be updated and the densities of the next timestep
pi(t+At) are calculated.
Since the method is Lagrangean by nature it is
straightforward to allow for nonNewtonian liquids with
thixotropy all computational nodes move around with
fluid velocity and can thus carry their own history. In the
experiments it was found that the sealing material studied
behaved like a Bingham fluid, that is, that the apparent
viscosity of the material could be expressed like:
TZi /yield
Y =t ,U (6)
at least in the stationary limit. The theory of this type of
material response to stress was first developed by Eugene
Bingham and a good description of the dynamics of plastic
fluids can be still be enjoyed in his textbook from 1922. For
the transient part of the dynamics it is assumed that the
material is of first order, obeying a simple ODE expressed
as:
ui /i i
at Tthi,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
impact of sealing droplets is so short in comparison to the
other timescales of the flow, this process could not be
resolved and it was therefore necessary to introduce an
impact coefficient of restitution between impacting SPH
material and wall boundaries.
The same inlet conditions, boundary conditions and material
properties are used for the simulations of the moving
applicator. In order to keep memory requirements
reasonable, SPH droplets were frozen after 1.0 s, which was
considered long enough for the flow to settle with an the
applicator moving at a velocity of 250 mm/s and cone
diameter at point of impact of maximum 40 mm, the direct
flow from the applicator itself will have passed any position
along its trajectory after less than 0.1 s, leaving 0.9 s for the
flow to settle before the droplet is frozen.
Results and Discussion
From the rotary rheometer experiments, yield stress,
Bingham plastic viscosity and structure modification
timescales were determined. Apparent viscosity as a
function of time is shown in Figure 6 and the stationary
shear stress is plotted versus rate of deformation in Figure 7.
Rheometric test series
500 1000 1500 2000 2500
Time [s]
3000 3500
Figure 6: Apparent viscosity as a function of time,
following the stepwise change in shear rate that are shown
(7) in Figure 3.
Simulations were set up to resemble all three experiments;
rotary rheometer, stationary applicator and moving
applicator. In order to numerically duplicate the rheology
experiments a flat cylinder of SPH particles was set up
between two walls. The boundary condition with respect to
the top wall was set to the rotational speed as specified in
the rheometer experiment; the lower wall boundary was
kept fixed. In order to compute the torque in the case of the
rotary rheometer an integration of the torques from the top
layer of the SPH particles were calculated.
In the stationary applicator simulation, SPH particles were
injected with velocity and direction towards the target plate
according to the corresponding measurements. Material
properties were taken from the numerical and physical
experiments of the rotary rheometer. As the timescale in the
300
200
100
0
0 100 200 300 400 500 600 700 800
Shear stress [Pa s]
Figure 7: Stationary shear stress versus rate of
deformation, rheometer experiments (red triangles) and
numerical results (blue circles).
U.~=Pi (t+ t =ui IVPt+ t
In Figure 8, the different transients as response to the
stepwise changes in settings are shown normalized in the
same plot. For the higher rates of deformation the difference
in viscosity was too small to reliably resolve in this analysis,
so only rates of deformation up to 50 s' could be used in
determining the structure adaptation timescale. The rates of
deformation in the actual material can be estimated using
impact velocity and string thickness as:
U
?_ = (8)
D
With characteristic tangential velocity U 1 m/s and sealing
film thickness D 1 mm the rate of deformation at the
points of impact can thus be as high as 1000 s1, but only for
very short periods of time.
* Shear rate Increased to 1 /s
* Increased to 5
A increased to 10
 increased to 30
* increased to 50
 Exponential cure fit
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
and to verify the physics of sealing material impacting the
plates.
A snapshot from footage of the stationary applicator
spraying high above the target plate is shown in Figure 9.
From measuring the length of the black streaks, which in
reality are sealing droplets, and from knowing the shutter
time on the camera the injection velocity of the droplets was
estimated to 10 m/s.
decreased to 50
decreased to 30
decreased to 10
decreased to 5
Shear rate decreased to 1 /s
1.2
1
0.8 A
0.4
Figure 8: Rheology transients. The average time scale for
0.2 .
be questioned. On the other hand, as the estimated
thixotropic timescale of around 40 seconds is rather long
(50 10 ,0 2(1
0.2
Time [s]
Figure 8: Rheology transients. The average time scale for
viscosity adaptation was estimated to 40 s.
As can be seen in Figure 8, the shape of the transients does
not conform very nicely to the first order exponential curve,
meaning that the first order assumption for the transient may
be questioned. On the other hand, as the estimated
thixotropic timescale of around 40 seconds is rather long
compared to the other timescales in the flow simulations
the assumption still has some merit. The changes in
viscosity during the passage of the applicator are
nevertheless large enough to have an impact on the
simulation result. In numerical experiments using a fixed
viscosity instead of the first order assumption it has not been
possible to reproduce features like walls and impact ditches
seen in both experiments and the more advanced
simulations.
A simulation of the rheometer experiment was first set up
and run to verify that the fundamental models of the
material theology did give the expected results. In the
following simulations, a hollowcone nozzle was
implemented. The first simulations using this nozzle were
made to duplicate the experiments with the stationary nozzle
Figure 9: Stationary applicator with sealing droplets
visible. Millimeter scale in the background.
In Figures 10 and 11, results from the stationary applicator
experiments and simulations are summarized. It can be seen
that the width of the applied strings follow the experiments
really well but that there is some discrepancy in the
comparisons of bead height. When visually comparing the
numerical and physical experiments, see Figures 3 and 12,
the shape of the beads are matching. As both width and
shape are similar there must be a discrepancy in bead
volume, and as the SPH simulations are inherently mass
conserving the only possible conclusion is that the volume
of applied material in the experiments somehow is larger
than in the simulations. As the mass of the applied material
is the same in experiments and simulations there seems to
be some entrainment of air as the sealing droplets impact the
target, which possible could make the material form a foam
upon impacting the target. After cutting one of the larger
beads open, small bubbles could be observed on the inside
which possibly could support this hypothesis. Further, the
hardened material is slightly elastic, which may also be an
indication of a small scale spongy structure being present. A
photo of the cut up material is shown in Figure 13.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0 08
007
006
005
E 004
n 003
002
001
0
u uuo 
e
0 0050
0 004
E 0003
 0 002
I
s .
4
0 . ..00 1. .
0 01 02 03 04 05 06 07
Spray time [s]
Figure 11: Bead height versus time, stationary applicator.
,z~: K'
A
Figure 12: Simulated sealing application, using the same
settings as in the case shown in Figure 3. Note the wall
buildup of material surrounding an impact ditch, similar to
what is seen in the physical experiments.
Figure 14: Comparison of experimental (left) and simulated
(right) sealing application. Top, middle and bottom images
show application of strings with applicator to plate distances
of 20, 40 and 60 millimeters respectively.
008
007
006
" 005
0 004
2 003
002
001
0
Figure 13: Photo of cut open sealing bead. Small cavities
are visible to the naked eye.
0 001 002 003 004 005 006 007
Nozzleplate distance [m]
Figure 15: Bead width versus applicator to plate distance.
A Exper ments e SPH In Figure 14, a rendition of a string application is shown
together with photos from the experimental sessions. The
_similarity to the physical experiment with the same settings
n is striking. To quantitatively compare simulations to
experiments, bead width as a function of applicator to target
distance was measured and compared. The results are
shown in Figure 15. As can be seen, for the shorter
applicator to target distances the results coincide really well
but in the case of larger distances the difference is
0 01 02 03 04 05 06 07 increasing. One explanation for this may be that airliquid
Spray time [s] interaction is completely disregarded in the simulations.
0: Bead width versus time, stationary applicator. This assumption holds true to a great extent as the transit
;rimental points at 0.3 and 0.5 seconds are partially time from applicator to target is small, but as the distance
behind the numerical results, increases there is more time in which drag forces can act on
the sealing material, reducing impact velocity and
1 subsequently reducing the horizontal flow of material on the
A Experiments SPH plate thus giving a narrower bead. As the measurements of
__droplet velocity were rather crude, this hypothesis could not
be verified using the measurement techniques available.
Figure 1
The expe
hidden b
0 008
0007
0006
A Experiments SPH ...... Linear fit
Sit
__ _ _ ___ *'
The simulation time for the typical simulations of about two
seconds of physical time was six hours using one processor
on a Intel Core2 Quad CPU with 2.40GHz clock frequency
and a time step of 106 s. The number of active SPH
particles in these simulations were at maximum 10 000.
Conclusions
Simulations and experiments of the process of applying
sealing materials on car bodies have been performed. The
results from the rotary rheometer experiments show that the
sealing material behaves like a Bingham fluid, at least in the
stationary limit. The results from the simulations show good
agreement with the experiments using the hollowcone
nozzle applicator. Although the work presented here only
has been concerned with the materials and applicators
specific to the sealing process in a car paint shop the results
and the methods are general and should apply to other (non
turbulent) flow scenarios with other viscoelastic fluids.
Acknowledgements
The authors wish to acknowledge the financial support by
the Swedish Governmental Agency for Innovation Systems,
VINNOVA. The authors are also grateful for the valuable
support from our industrial collaborating partners Volvo Car
Corporation and Saab Automobile AB. Photos in Figure 1, 3,
4, 9 and 14 by courtesy of Volvo Car Corporation.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
References
E.C. Bingham (1922), Fluidity and Plasticity, McGrawHill,
New York, 1922
S.J. Cummins and M. Rudman (1999), An sph projection
method, J. Comp. Phys, vol 152, 584607
M. Ellero and R.I. Tanner (2005), SPH simulations of
transient viscoelastic flows at low Reynolds number, J.
NonNewtonian Fluid Mech, 132, 6172
J. Fang, A. Parriaux, M. Rentschler and C. Ancey (2009),
Improved SPH methods for simulating free surface flows of
viscous fluids, Applied Numerical Mathematics, Vol 59, 2,
251271
J.J. Monaghan (2005), Smoothed particle hydrodynamics,
Rep. Prog. Phys, 68, 17031759
A. Rafiee, M.T. Manzari and M. Hosseini (2007), An
incompressible SPH method for simulation of unsteady
viscoelastic freesurface flows, International Journal of
NonLinear Mechanics, 42, 12101223
S. Shao and E.Y.M. Lo (2003), Incompressible SPH method
for simulating Newtonian and nonNewtonian flows with a
free surface. Advances in Water Resources, 26, 787800
