A Fluid Resistance Map Method for Real-time Haptic
Interaction with Fluids
Yoshinori Dobashi

Shoichi Hasegawa

Mitsuaki Kato

Hokkaido University
Kita-ku, Kita 14, Nishi 9
Sapporo 060-0814, Japan+81.11.706.6530

Tokyo Institute of Technology
4259 Nagatuta
Yokohama 226-8503, Japan+81.45.924.5050

The University of Tokyo
7-3-1, Hongo, Bunkyo-ku
Tokyo 113-8656, Japan+81.4.7136.3942

doba@nis-ei.eng.hokudai.ac.jp

hasegawa@pi.titech.ac.jp

kato@fel.t.u-tokyo.ac.jp

Makoto Sato

Tsuyoshi Yamamoto

Tomoyuki Nishita

Tokyo Institute of Technology
4259 Nagatuta
Yokohama 226-8503, Japan+81.45.924.5050

Hokkaido University
Kita-ku, Kita 14, Nishi 9
Sapporo 060-0814, Japan+81.11.706.6530

The University of Tokyo
5-1-5, Kashiwanoha
Kashiwa 227-8561, Japan+81.4.7136.3942

msato@pi.titech.ac.jp

yamamoto@ist.hokudai.ac.jp

nis@is.s.u-tokyo.ac.jp

ABSTRACT

1. INTRODUCTION

Haptic interfaces enable us to interact with a virtual world using
our sense of touch. This paper presents a method for realizing
haptic interaction with water. Our method displays forces acting
on rigid objects due to water with a high frame rate (500 Hz). To
achieve this, we present a fast method for simulating the
dynamics of water. We decompose the dynamics into two parts.
One is a linear flow expressed by a wave equation used tocompute water waves. The other is a more complex and nonlinear flow around the object. The fluid forces due to the nonlinear flow is precomputed by solving Navier-Stokes equations,and stored in a database, named the Fluid Resistance Map. The
precomputed non-linear flow and the linear flow are combined tocompute the forces due to water.

Intuitive and natural interactions with virtual world are key
mechanisms for enhanced perception of the virtual worlds. The
haptic rendering is one of such technologies and many methods
have been presented [16] [17] [18]. In a typical virtual reality
system with a haptic interface, the user controls the movement ofa virtual object using the haptic interface, and the systemcomputes forces when the object contacts with other objects. The
user can also feel the resulting forces via the haptic interface. This
provides the user with a natural interactions with the virtual world.
Previous methods have efficiently displayed forces between rigid
objects and forces due to soft objects [3] [11] [18] [19]. In thispaper, we focus on haptic interactions and rendering with fluids,
especially water. In computer graphics, while several publications
exist on the visual simulation of fluids such as smoke, water [5][9] [10], there are few published techniques for the haptic
interaction and rendering with fluids. Examples of the forces due
to interaction with fluids are aerodynamic resistances on an
airplane, hydrodynamic resistances on a ship, water buoyancy,
frictional and pressure resistances with the oars of a boat, and alure used in fishing.

CR Categories
I.3.7 [Computer Graphics]: Three-Dimensional Graphics andRealism - Animation; Virtual Reality; I.6.5 [Simulation andModeling]: Model Development – Modeling methodologies; I.6.3[Simulation and Modeling]: Types of Simulation - Animation.

General Terms
Algorithms.

Keywords
Virtual Reality, Haptics, Fluid
Computational Fluid Dynamics.

Resistance,

Simulation,-----------------------------------------------------Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that
copies bear this notice and the full citation on the first page. To copy
otherwise, or republish, to post on servers or to redistribute to lists,requires prior specific permission and/or a fee.
VRST'06, November 1–3, 2006, Limassol, Cyprus.
Copyright 2006 ACM 1-59593-321-2/06/0011...$5.00.

In this paper, we present a method for very rapidly estimating anddisplaying forces acting on a virtual object due to water. In thefield of computational fluid dynamics and ship design engineering,
researchers have developed methods to numerically estimate theresistance due to fluids by solving the Navier-Stokes equations [1][6]. However, their estimation techniques while numericallyaccurate are prohibitively time consuming. This is especially truefor haptic rendering, where a high frame rate is required.
We additionally present a water simulation model for haptic
interactions. We decompose dynamics of water into two parts, alinear flow and a non-linear flow. The linear flow is the flow far
from the object. The linear flow is well approximated by wave
equations derived by linearizing Navier-Stokes equations under
the assumption that the flow of water is slow. However, thisassumption is not valid for the flow close to the object. The flow
around the object is complex and highly non-linear. This nonlinear flow has to be computed by directly solving the NavierStokes equations, which unfortunately is compute intensive. So
we precompute the non-linear flow around the object and create a
database, which we call a fluid resistance map, or FRM for short.

Figure 1: An example application of our method, a virtual canoe simulator. The user can row the virtual canoe using the haptic
this system. The forces and the torques on the paddle are calculated in real-time by our method. The user feels the rowing motion
via the haptic interface. The motion of water is also calculated by our method and displayed in real-time.
In computing FRM, we separate the effect of the gravity, which
creates buoyancy, from the total water force. This eliminates thedependency of FRM on the depth from water surface.
Consequently, the presented method computes three types offorces: 1) the buoyancy, 2) the frictional resistance due to thefriction between the object surfaces and water, and 3) the pressure
resistance due to the movement of the object. In the following, wecall the cumulative force of the frictional and the pressure
resistances, as fluid resistance. FRM stores the fluid resistanceacting on points sampled on the object surface.
Fig. 1 shows an example application of our method, a virtual
canoe simulator. The user can control the paddle of the canoe, andfeel the forces and torques calculated by our method at 500 Hz.
As shown in this example, when we try to move the object in thewater, we often feel not only the forces against the movement
direction but also the forces perpendicular to the movement
direction. This is due to the non-linear flow around the object.
That is, vortices generated behind the object cause unexpectedfluctuations of the fluid resistance. The fluid resistances are timevarying because of the generation of these vortices. To simulate
this, FRM stores the time-varying fluid resistance.
Our method has the following features:(1) Decomposition of the water dynamics into linear flow andnon-linear flow for real-time force display.
(2) Computation and force feedback of the realistic force taking
into account the buoyancy, frictional resistance, andpressure resistance.
(3) Simulation of the time-varying fluid resistance.
(4) Rendering of force field due to water by coupling the linear
flow (wave equation) and the non-linear flow (FRM).
We demonstrate our method by using a 6 DOF (degrees-offreedom) haptic interface. The user perceives the torques as wellas the forces.
The rest of the paper is organized as follows. First, in Section 2,
we discuss the previous work. Next, in Section 3, we present oursimulation method for the dynamics of water. In Section 4, we
provide an overview of the implementation of the simulation
method. In Sections 5 and 6, we provide our method for creating
an FRM and the method for the real-time display of the forces. In
Section 7, we verify our method by comparing simulation withmeasurement made with our experimental setup. Furthermore,
several examples are demonstrated using haptic interfaces. Next,

Section 8 is devoted to discussion of the use of our method.
Finally, Section 9 concludes this paper and discusses future work.

2. PREVIOUS WORK
There have been several published methods for haptic interaction
and rendering. Mark et al. addressed issues in adding force
feedback to graphics systems [16]. Ruspini et al. presented ahaptic interface library “HL” [20]. McNeely et al. presented amethod for the haptic rendering of complex virtual environments
by voxelizing objects [17]. Otaduy and Lin presented asimplification algorithm for faster collision queries with sufficientaccuracy for the haptic rendering [18]. Hasegawa and Satopresented a method for real-time rigid body simulations for haptic
interactions [11]. These methods realize realistic haptic
interactions with virtual worlds. These methods, however, cannot
handle forces by objects interacting with fluids. Durbeck et al.
displayed forces due to interactions with fluid based on the resultof computational fluid analysis [8]. However, they simply
displayed the forces according to the flow field and the forceswere not physically-based.
To our knowledge, Baxter and Lin tried for the first time todisplay the forces due to fluids based on physical phenomenon [2].
They succeeded in a real-time haptic display of fluids in twodimensions by directly solving Navier-Stokes equations. However,computing the fluid force in three-dimension by solving NavierStokes equations is time-consuming and therefore the real-time
display of the force becomes difficult. In addition, an object
controlled by the user is discretized by using grids. The
discretized shape changes from frame to frame when the object
moves. This causes undesirable noise in the resulting force.
In the field of computational fluid dynamics (CFD), ship design
engineering and aeronautics, several researchers have developed
methods for computing aerodynamic or hydrodynamic resistances[1] [6]. However, their purposes are to obtain optimal design of an
airplane or a ship, that minimizes the drag resistances. Thesemethods require very expensive computational calculations.
Therefore, they are not suitable for real-time applications such as
interactions with virtual worlds.

3. OUR SIMULATION METHOD OF
WATER DYNAMICSFirst, we present a method for simulating water dynamics for
haptic interactions. Next, we describe a model for thecomputation of forces between an object and surrounding water.

3.1 Water Dynamics for Haptic Interactions
An object in the water receives forces due to water pressure andfrictional drag. We assume that pressure can be decomposed into
two components, that is, the pressure due to the movement of theobject and the pressure due to gravity. The pressure due to thegravity causes buoyancy. As a result, the presented model treats
three types of forces: 1) buoyancy, 2) frictional resistance, and 3)
pressure resistance. To display these forces in real-time, we
decompose the dynamics of water into linear and non-linear flows
as shown in Fig. 2(a). The linear flow is used to simulate thewater flow far from the object. The non-linear flow is used forcomplex flow near the object. The water flow far from the object
is generally slow. In this case, the water dynamics is wellapproximated by linearizing Navier-Stokes equations, resulting in
two-dimensional wave equations [13]. On the other hand, the nonlinear flow around the object is complex and therefore simulated
by solving Navier-Stokes equations. The linear flow can beobtained in real-time since the computational cost for the wave
equation is small. However, solving Navier-Stokes equations inreal-time to obtain the non-linear flow is difficult. Therefore, nonlinear flow is precomputed for various fluid flows and stored as
an FRM.
The linear flow (wave equations) and the precomputed non-linear
flow (FRM) are combined for real-time haptic interactions. The
idea is shown in Fig. 2 (b). First, the shape and velocity of water
waves are obtained by solving the wave equations. The shape of
the water wave is used to compute the buoyancy. The velocity isused to compute the relative velocity of the flow against thevelocity of the object, which is controlled by the user. The fluid
resistance is then calculated using the relative velocity byreferring to the FRM. The sum of the buoyancy and the fluid
resistance is displayed to the user. In addition, the pressure due tothe non-linear flow stored in the FRM is applied to the wave
equations as external forces. This generates new waves. In thisway, the linear flow and the non-linear flow are combined todisplay the force and simulate the dynamics of water efficiently.

3.2 Force ModelIn this subsection, our computation model for the buoyancy andthe fluid resistance are described by using Fig. 2(c).
The buoyancy is calculated approximately by using hydrostatic
water pressure. The hydrostatic water pressure is the pressure
under static water. The force due to hydrostatic pressure, fb(Q), on
point Q acts in the direction of the normal of point Q (see Fig.
2(c)) and its magnitude is proportional to the depth of point Q
below the water surface. That is,fb (Q)

 U w wQ n Q ,linear flow

objectvortexwater surface
non-linear flow(a) Linear flow and non-linear flow.
buoyancy linear flow
(wave equations)watervelocity+

The frictional and the pressure resistances are caused due to theviscosity of water when the object is placed in a flow or when theobject moves under water. Basically, these resistances act toprevent the movement of the object. In the following, we assume
that the relative velocity of the flow to the velocity of the object isu (see Fig. 2(c)).

relative
velocity non-linear flow(precomputed FRM)
objectvelocity

haptic interface

fluid resistance(b) Block diagram for our simulation.
object motion
wQ

nQ

frictionz*Qrhydrostaticpressurepressurewater surface
objectOrg
center ofmass
gravityvortexrelative velocity of flow u(c) Force model.
Figure 2: Our simulation method for water dynamics.

The frictional resistance is due to the friction between the object
surface and water. To compute the frictional resistance accurately,
we have to compute the flow in the so-called boundary layer,
which is turbulent flow layer very close to the object surface [1][12]. However, the flow in the boundary layer is very complicatedand numerical analysis is difficult. Therefore, we employ one of
the simplest approximations. That is, the frictional resistance ff(Q)on point Q is proportional to the gradient of the relative velocity
with respect to the normal direction of point Q [6].

(1)

where Uw is the density of water, wQ is the depth of point Q, nQ isthe normal vector of point Q. The buoyancy force of the object isobtained by integrating the above equation over the submerged
part of the surface of the object. By definition, the hydrostaticpressure is the water pressure when there are no waves. However,
we use depth wQ taking into account the waves as shown in Fig.
2(c). This allows us to simulate fluctuations of the buoyancy due
to the waves.

-pressuref f (Q)

P

wu*
wz *

,

(2)

z* 0

where P is a dynamic viscosity of water and u* is the velocity in a
local coordinate system whose z* axis is the normal vector ofpoint Q (see Fig. 2(c)). The frictional resistance is small for fluidswith small viscosity. However, it becomes important for fluidswith large viscosity.
The pressure resistance is caused by the difference between
pressure in front of the object and pressure behind the object. In
general, the pressure behind the object is lower than that in frontof the object. The pressure resistance fp(Q) on point Q isexpressed by the following equation.
f p (Q ) pn Q ,(3)

flow ene4

haptic source

flow direction
e3

haptic region

1

2

4

5

flow e2
flow e1e2

3
6

relative velocity
haptic source
pointer object
water
surfacee1pressureRp(l, e1, s)

pointer object

frictionRf(l, e1, s)(b) Fluid resistance map for haptic source l.

(a) Numerical fluid analysis (preprocess).

hydrostaticpressure

fluid resistance(c) Real-time process.

Figure 3: Overview of implementation.
where p is the pressure on point Q.
The total force acting on the object is obtained by integrating thehydrostatic pressure fb (Eq. 1) and the fluid resistance ff and fp(Eqs. 2 and 3) over the underwater part of the object surface, andthen adding the effect of gravity of the object. That is,F³ (f (Q)  f (Q)  f (Q))ds  mg ,
S

b

p

(4)

f

where S indicates the underwater part of the surface, m is the massof the object, and g is the downward acceleration vector of gravity.
A torque is also generated according to the above forces. The
torque at point O in the object is calculated by:T³ {r u (f (Q)  f (Q)  f (Q))}ds  m(r u g) ,
S

b

p

f

g

(5)

where r is a vector from point O to point Q, rg is a vector frompoint O to the center of mass (see Fig. 2(c)).
For computing the fluid resistance (ff and fp), the pressure and thevelocity distribution around the object are required. Moreover, thefluid resistance is time-varying due to the generation of vortices
behind the object, such as the Karman vortex (see Fig. 2(c)). To
simulate this efficiently, we make use of the following
experimental properties for fluid resistance [4] [7] [15].
1. The magnitude of the pressure resistance is proportional to
the square of the relative speed of the flow.
2. The magnitude of the frictional resistance is proportional to
the relative speed of the flow.
3. The frequency of the fluctuation of the fluid resistance (i.e.
the frequency of the generation of the vortices) is
proportional to the relative speed of the flow.
Once the fluid resistance at base speed u0 is obtained, we cancompute the fluid resistance of an arbitrary speed by using theabove three properties.

4. OVERVIEW OF IMPLEMENTATION
Fig. 3 shows the overview of the implementation of oursimulation method presented in the previous section. First, the
FRM is created in a preprocessing step. Then forces and torques
are computed in the haptic rendering process. Objects other thanwater are assumed to be rigid and modeled by polyhedra. The
user controls a rigid object by using the 6 DOF haptic interface. In
the following, we call the user-controlled object the pointer object.
In the preprocess, the pressure and the velocity distributions
around the pointer object placed in various directions ei of
uniform flows are calculated by CFD as shown in Fig. 3(a). The

speed of the uniform flows is set to the base speed u0. Then, thefluid resistances are calculated by using Eqs. 2 and 3, respectively.
The results are stored in an FRM. As we mentioned before, thefluid resistance is time-varying. However, storing the timevarying fluid resistance for all points on the object surface
requires a large amount of memories. Fortunately, the fluid
resistances on neighboring points are very similar. We make useof this coherency. That is, as shown in Fig. 3(a), the surface of thepointer object is subdivided into several regions. We call theseregions haptic regions. At the center of each haptic region, we
generate a point haptic source. FRM stores the fluid resistanceacting on the haptic region. Fig. 3(b) shows contents of FRM forpoint haptic source l. FRM consists of the pressure resistance Rp(l,ei, s) and the frictional resistance Rf(l, ei, s). l is the haptic source
number and s indicates time. Time series of xyz components ofthese resistances are stored for each direction of the flow ei.
In the haptic rendering process, the buoyancy and the fluid
resistance are calculated in real-time. The water surface is also
displayed in real-time by solving the two-dimensional wave
equations. The water surface is expressed as a height field
represented by a two-dimensional mesh. At each grid point of themesh, the height and the velocity of the waves are calculated.
Then, the underwater portion of the object surface is extracted.
The hydrostatic pressure at each point on the underwater part iscalculated to obtain the buoyancy. The fluid resistance iscalculated by referring to FRM using the relative velocity of thepointer object against the wave velocity (see. Fig. 3(c)). The
obtained fluid resistance is used not only for the force display but
also for generating waves due to the movement of the object.
Details of the preprocessing step and the haptic rendering process
are described in the subsequent sections.

5. FLUID RESISTANCE MAP
In this section, we describe the method for computing FRM. FRM
is created by the numerical analysis of the flow around the pointerobject placed in uniform flows with various directions ei. The
speed of the flows is set to the base speed u0.

5.1 Numerical Fluid Analysis
As shown in Fig. 4(a), an analysis region is specified around thepointer object and is divided into voxels. The flow inside thisanalysis region is simulated by a finite difference method [21].
The number of voxels is determined experimentally so that theshape of the object is well approximated by the voxels. The time
step for the fluid simulation is set as 1/500 sec. Velocity u = (ux,uy, uz) and pressure p are assigned to each voxel. As a boundary
condition, the velocity at y = 0 is (0, u0, 0). A natural boundaryseries into the frequency domain by FFT and removes highfrequency components, whose amplitudes are smaller than aspecified threshold. After that, we transform them back into thetime domain. This process makes FRM periodic.

pointer object

6. REAL-TIME HAPTIC RENDERING
flow(0, u0, 0)zyx(a) Fluid analysis.
frictionIn the real-time process, the images of water and the objects are
displayed at 30 Hz. On the other hand, the forces are computed at
500 Hz. Experimentally, we found that 500 Hz is sufficient for
stable calculations since the forces due to water do not change
rapidly. The following procedure is executed at every 1/500 sec.
1. Update of the position and the velocity of the pointer object
using the pointer position of the haptic interface.

pressure

2. Computation of the buoyancy
3. Computation of the fluid resistanceQ

4. Computation of the water surface
5. Display of the resulting forces and the torques using thehaptic interface

object surface(b) Computation of pressure and friction.
Figure 4: Creating fluid resistance map.
condition is applied to other boundaries. That is, the differentials
of velocity and pressure normal to the boundaries are 0. The
initial velocities and pressures of all voxels are zero.
Navier-Stokes equations are solved numerically under the above
conditions. Then the frictional and the pressure resistances arecalculated by using Eqs. 2 and 3. In this fluid analysis, thebuoyancy is not taken into account since it is calculated in thereal-time process using Eq. 1. This is equivalent to solving
Navier-Stokes equations without taking into account the gravity.
The fluid simulation is computed for a specified number of time
steps. During the early time steps, however, the simulation is in atransition state and is unstable. So, the computation of the fluid
resistance begins after the simulation becomes stable.

5.2 Generating Fluid Resistance Map

Before starting the simulation, sample points to detect theunderwater part of the pointer object are generated on the surfaceof the pointer object as follows. Each polygon of the pointerobject is subdivided into a specified number of triangles. Then thesample points are generated at the centers of the triangles. An area
and a normal of the triangle are stored for each sample point.
In the following subsections, Details of steps 3, 4, and 5 in theabove procedure are described.

6.1 Computing Buoyancy
The buoyancy is obtained by computing the hydrostatic pressures(Eq. 1) at the sample points on the surface of the pointer object.
First, the sample points under the water surface are extracted.
Since the water surface is represented by the height field, it isvery easy to check if the sample point is under the water surface
or not. Then the hydrostatic pressures are calculated at thosesample points. The buoyancy is obtained by summing theresulting hydrostatic pressures after multiplying the areas stored
for the sample points.

6.2 Computing Fluid Resistance

The frictional and the pressure resistances are calculated by using
velocities and pressures at voxels corresponding to the object
surface. As shown in Fig. 4(b), the object is represented by voxels(Fig. 4(b) shows a two-dimensional case for simplicity).
Therefore, the resistances on the voxel are obtained by computing
the resistances acting on the side faces of the voxel. The pressure
resistance is perpendicular to the faces and the frictional
resistance is tangential to the faces as shown in Fig. 4(b). Then,the sum of the resistances of the voxels corresponding to thesurface of each haptic region is stored.

The fluid resistance is calculated by using FRM. FRM stores thefrictional resistance Rf(l, ei, s) and the pressure resistance Rp(l, ei,s) acting on haptic source l of the pointer object placed in a flow
(speed u0 and direction ei) at time s. To avoid confusion, we
denote the time in FRM as s. The actual time is t. The fluid
resistance acting on the pointer object moving at an arbitrary
velocity is calculated by making use of the three properties
described in Section 3.2. The basic idea is to interpolate theresistances stored in FRM by using the velocity of each haptic
source. Details are described in the following.

The above computation is repeated at each time step of thesimulation. After the simulation, time series of the fluid
resistances are stored in FRM (see Fig. 3(b)). Unlike the methodby Baxter and Lin [2], our method does not produce theundesirable noises. The fluid resistances stored in the FRM aresmooth and continuous since the pointer object is discretized onlyin this preprocessing step and the descretized shape does notchange throughout the simulation. In the real-time haptic
rendering process, we use FRM periodically since the size of
FRM is finite. However, FRM obtained by the above method is
not periodic and hence the user would feel the discontinuous
forces at the ends of FRM. To address this, we transform the time

Let us assume that time t = 0 when the haptic rendering process
starts. The fluid resistance is calculated at a short time interval 't(= 1/500 in our case). The total fluid resistance is obtained bysumming the resistances at each haptic source. In the following, tk= k't, where k is a non-negative integer indicating time steps inthe simulation. We explain the computation method for haptic
source l.
First, a ratio of area of the underwater part of the haptic region of
haptic source l is calculated by using the sample points generated
on the surface of the pointer object. Let us denote the ratio as Dl.
When Dl = 0, the fluid resistance is set to 0. Otherwise, first, thevelocity of haptic source l is calculated according to theobject. In addition, the object pushes the water with pressure p
due to the movement of the object. This effect is simulated asfollows. Let us consider a pressure at point Q(Qx, Qy, Qz) on the
pointer object as shown in Fig. 5. The total pressure pw at point Q
is the sum of pressure p due to the movement of the pointer object
and the weight of water above point Q. That is,water surface z=h(x, y)

Uwg(h-Qz)ppw = Uwg(z-Qz)-p
Q(Qx,Qy,Qz)

pw

pointer object

wu
wt

movement direction
Figure 5: Computing water surface.
movement of the pointer object. The velocity of the waves of thewater surface is also calculated by the method described in thenext subsection. The relative velocity of the haptic source against
the velocity of the water waves is used to calculate the fluid
resistance. This can simulate the reaction due to water waves. Let
us denote the magnitude of the relative velocity as ul(tk) and the
direction as el(tk). Then the fluid resistance is calculated byreferring to FRM using ul(tk) and el(tk). That is, the pressure
resistance fl,p and the frictional resistance fl,f are calculated by thefollowing equations so that the three properties in Section 3.2 aresatisfied.

­f l , p (t k ) D l (ul (t k ) / u 0 ) 2 R p (l , e l , sk )°°(6)
®f l , f (t k ) D l (ul (t k ) / u 0 ) R f (l , e l , s k ) ,°°¯sk (ul (t k ) / u 0 ) 't  sk 1
where s0 = 0. As shown in the above equations, when the relative
speed is ul(tk), the value at intervals of (ul(tk)/u0 't in FRM is used
to compute fl,p and fl,f. This is equivalent to making the frequencyof the fluctuation of the resistances (ul(tk)/u0) times higher (orlower). The magnitudes of the pressure and the frictional
resistances are proportional to (ul(tk)/u0)2 and (ul(tk)/u0),
respectively. So, the properties described in Section 3.2 aresatisfied.
The above processes are repeated for all haptic sources and thetotal fluid resistance is obtained by summing them.

6.3 Computing Water Surface
When the pointer object moves near the water surface, waves are
generated. We simulate this by combining FRM with the method
proposed in [13].
Similarly to [13], we employ the following assumptions: theeffect of the viscosity is very small, the speed of the water waves
is slow, and amplitudes of the waves are very small to the depth
from the bottom of the water. Then, the water surface isrepresented by a height field h(x, y). Navier-Stokes equations arelinearized as follows.
wu
wt



1

Uw

p w ,

U w g (h  Qz )  p ,(9)where g is the acceleration of gravity. By putting Eq. 9 into Eq. 7,
the following equation is obtained.

(7)wh d  u ,(8)wt
where u = (ux, uy, uz) is a velocity vector, Uw is a fluid density
coefficient, d is a parameter of the wave equation, meaning the
average water height from the bottom. The waves caused by the
motion of the pointer object are simulated as follows. Inside thepointer object, the velocity of water is forced to the velocity of the gh 

1

Uw

p .

(10)

We solve Eqs. 8 and 10 to obtain the shape and the velocity of thewater waves. These equations are solved numerically by a finite
difference method. To compute the water surface efficiently, we
assume that the vertical speed of the waves is very small
compared to the horizontal speed, that is uz = 0. The pressure p is
obtained by a weighted average of pressures at neighboring haptic
sources, calculated by the method described in Section 6.2. We
use a Gaussian function for computing the weights. We consider
only the upper part of the object surface as shown in Fig. 5. The
pressure on the upper part of the object surface directly reaches
the water surface under the assumption of small viscosity. Since
the water surface is represented by two-dimensional grids, weevaluate the pressure gradient in Eq. (10) at each grid point above
the surface of the pointer object.

7. RESULTS
In this section, we first verify the validity of our method bycomparing measured data with simulation results. Next, several
examples are demonstrated, where the 6 DOF haptic interface
named SPIDAR-G [14] is used to interact with the virtual world
while display the forces and the torques. We connected SPIDARG with a desktop PC (Pentium IV 2.8GHz).

7.1 Comparison with Measured Data
We measured the forces acting on a parallelopiped. The
dimension of the parallelopiped is 5 cm x 1 cm x 50 cm. We
chose the parallelopiped in order to verify the capability of our
method to handle the anisotropic shape. The experimental setup is
shown in Fig. 6. We used a robot arm to measure the force. A
force sensor is attached to the tip of the arm. The parallelopiped isthen attached to the force sensor. We programmed the arm to
track a quarter of a circle with the 4/5 part of the parallelopiped
underwater. We chose two speeds of the arm movement, 7.5 cm/s,
and 10 cm/s. Furthermore, for each speed, we measured the forces
twice by attaching the parallelopiped at different angles, that is, 0
and 60 degrees. Zero degree corresponds to moving theparallelopiped in the direction perpendicular to its wider face. We
measured the forces in the inverse and perpendicular directions to
the movement of the arm. We compared the measured data withthe simulated forces under the same situation. For the simulation,
we created FRM by two-dimensional fluid analysis using thecross-sectional shapes of the parallelopiped. The number of thedirections of uniform flows for creating FRM is 36. The basespeed is u0 = 50 cm/s.
Fig. 7 shows our comparison. Figs. (a) and (b) correspond to theforces at speed 7.5 cm/s, respectively and Figs. (c) and (d) to theforces at speed 10 cm/s. We measured the forces by making therobot arm go back and forth several times. The sharp peaks in Fig.
7 correspond to the turning point of the movement. The blacklines correspond to the measured forces and the gray lines to thesimulated forces. The thick lines correspond to the forces in themeasured30

15perpendicular direction0

simulated

perpendicular direction

0

-15

-15

measured30

force [N]15

force [N]robot arm

simulatedmovement direction

movement direction

force sensor

-30

-30
4.0

8.0

12.0

4.0
8.0
12.0
16.0 time [s](b) speed 7.5 cm/s, angle 60 degrees

16.0 time [s](a) speed 7.5 cm/s, angle 0 degrees
measured30

force [N]

15

15

0

0

-15

movement direction

-30
4.0

8.0

12.0

movement direction-30
16.0 time [s](c) speed 10 cm/s, angle 0 degrees
Figure 6: Experimental setup.

simulated

perpendicular direction

perpendicular direction-15

parallelepiped

measured30

force [N]water

simulated

4.0
8.0
12.0
16.0 time [s](d) speed 10 cm/s, angle 60 degrees

Figure 7: Comparison with measured forces.

direction of the movement and thin lines to the forces
perpendicular to the movement. At the turning point, there aresudden changes of the forces since the relative velocity of theparallelopiped to the surrounding water increases at that time.
Similar effects are also observed in the simulated data. As shownin Fig. 7, our method reproduced the forces very similar to themeasured forces. However, some differences are observed aroundthe peaks. This is because the inertia forces were included in themeasured data. When the direction of the movement changes, theforce sensor measured the forces from inertia of theparallelopiped as well as the forces due to water. The inertia
forces were not calculated in the simulation.

extending the SPIDAR-G haptic interface. The user controls thepaddle using this interface. The buoyancy and the fluid resistanceare computed not only for the paddle but also for the canoe. The
canoe moves according to the forces due to water as well as theforces due to the paddle. Fig. 10 shows a lure-fishing simulator.
To improve the reality, we attached to SPIDAR-G a stick thatsimulates a rod as shown in Fig. 10(a). The force on theunderwater fish lure is calculated and this drags the tip of the rod.
In these examples, the numbers of the haptic sources are 40 and
30 for Figs. 1 and 10, respectively. The water surface isrepresented by 200 x 200 and 160 x 200 grids for Figs. 1 and 10,
respectively.

7.2 Examples using Haptic Interface

In the examples shown in this subsection, the precomputation
times for FRM ranged from 1 to 24 hours. The sizes of FRM
ranged from 1 MB to 10 MB. In the haptic rendering process, theforces and torques are calculated and displayed at 500 Hz. Pleaserefer to the movies corresponding to each of the examples shownin this subsection, provided in the supplementary submitted files.

Figs. 8 and 9 show simple applications with the haptic interface.
In Fig. 8, the user controls a sphere with radius 5 cm and feel theforces in real-time when the sphere interacts with water. The
water surface is represented by 128 x 128 grids. The numbers of
haptic sources are 26 and 20 for Figs. 8 and 9, respectively. The
user moves the sphere in a clockwise circular orbit. Figs. 8(a)through 8(c) correspond to the images at time 0.7, 1.3, and 1.9
seconds. The arrows in these images indicate the direction of theforces. Fig. 8(d) shows the relation between the speed of thesphere and the force. Three components of the force are shown:forces in the movement direction, the horizontal direction, andvertical direction to the movement. The force against the
movement is dominant and the user hardly feels the forces
horizontal and vertical to the movement. This is because the flow
around the sphere is smooth and symmetrical. In Fig. 9, the user
controls a parallelopiped instead of the sphere. The dimension of
the parallelopiped is 8 cm x 1 cm x 30 cm. When the user tries tomove the parallelopiped forward as shown in this figure, thestrong fluid resistance acts in the horizontal direction (see Fig.
9(d)) they cause the parallelepiped to move left and rightalternatively. We often observe this phenomenon in the real world.
In Figs. 8 and 9, the images of the water surface are rendered by
mapping a sky image using an environmental texture mapping
technique.
Figs. 1 and 10 show more practical applications. Fig. 1 is a canoe
simulator. We constructed a human-scale haptic interface by

8. DISCUSSION
When using our method, the following points should be kept inmind.
First, our method cannot take into account randomness due toturbulent flow. However, we believe that our method reproduces
the distinctive features of the flow forces. Hence, our presented
method is very viable for haptic interfaces in a virtual reality
systems. To reproduce the randomness caused by turbulent flow,
Navier-Stokes equations must be fully solved taking into accountthe water velocity around the submerged object. However, this isvery time-consuming and therefore is not applicable to haptic
interactions requiring the high frame rate.
Secondly, when the user rotates the pointer object, the forcescomputed by our method are not totally correct. This is because
the FRM is created only for uniform flows. The forces due torotation are not stored in the FRM. In our method, the forces due
to rotation are calculated by assuming that each haptic source isplaced in a uniform flow at the speed of the rotation during theshort time interval of a real-time simulation. To simulate theforces due to the rotation accurately, the rotational motion of thepointer object must be simulated in a preprocessing step and theresulting forces need to be stored. Again, they will lead to asignificant increase in the computation time and the memory
requirement.
Thirdly, the proposed method cannot reproduce the forces due tothe other objects near the pointer object accurately. This effect issimulated approximately in the following way. The objects near
the pointer object can generate water waves calculated by wave
equations. These water waves affect the forces acting on the
pointer object since the force is determined by the relative
velocity of the pointer object against the velocity of the water
wave.
Finally, we would like to mention that our method can be applied
to the simulation of forces due to other fluids such as air. In thiscase, the buoyancy forces can be ignored since they are rather
small compared to the fluid resistance forces.

9. CONCLUSION
We have presented a method for the real-time display of theforces acting on rigid objects due to water. Our method candisplay the buoyancy, the frictional resistance, and the pressure
resistance. To compute these forces in real-time, we havepresented a method for simulating the dynamics of water. The
water dynamics was decomposed into linear flow and non-linear
flow. For linear flow, we used the wave equations to create water
waves. For non-linear flow, we have presented an FRM, the fluid
resistance map. The FRM is created in a preprocessing step bysolving Navier-Stokes equations, and the resulting fluid resistanceis stored. Then, the wave equations and the FRM are combined tocompute the forces due to water in real-time. The accuracy of our
method has been demonstrated by comparison with actual
measured data. The usefulness of the method has also been
demonstrated by several examples using the SPIDAR-G haptic
interface.
Our method provides a natural way to interact with fluids in a
virtual world. We believe that this greatly enhances the reality ofthe virtual environment.
There are a few things to be done in the future. Currently, we
determined the number of the haptic sources experimentally.
However, the computational cost for the real-time haptic
rendering process is proportional to the number. Therefore, an
adaptive method for determining the optimal number andpositions of the haptic sources need to be developed to speed up
the computation without losing the accuracy. Next, as we
discussed in Section 8, the fluid resistances corresponding tovarious velocities and rotational motions of the pointer objects
should be stored in FRM in order to display the forces accurately
under more complex situations. In this case, we have to develop acompression method suitable for FRM since the memory
requirement for storing such fluid resistances increases
significantly.

10. REFERENCES[1] Anderson, J. D., JR., Computational Fluid Dynamics,
McGraw-Hill (USA), 1995.
[2] Baxter, W. and Lin, M. C., Haptic Interaction with Fluid
Media, In Proceedings of Graphics Interface2004, Annual
Conference Series, 2004, 81-88.
[3] Baxter, B., Scheib, V., Lin, M. C., and Manocha, D., DAB,Interactive Haptic Painting with 3D Virtual Brushes, In
Proceedings of ACM SIGGRAPH 2001, Annual Conference
Series, 2001, 461-468.

[4] BloomerJ., Practical Fluid Mechanics for Engineering
Applications,” Marcel Dekker, 2000.
[5] Carlson, M., Mucha, P. J., and Turk, G., Rigid Fluid:
Animating the Interplay Between Rigid Bodies and Fluid,
ACM Transaction on Graphics (Proceedings of ACM
SIGGRAPH 2004), 23, 3, 2004, 377-384.
[6] Chung, T. J., Computational Fluid Dynamics, Cambridge
University Press, 2002.
[7] Dobashi, Y., Yamamoto, T., and Nishita, T., Real-timeRendering of Aerodynamic Sound Using Sound Textures
based on Computational Fluid Dynamics, ACM Transaction
on Graphics (Proceedings of SIGGRAPH 2003), 23, 3, 2003,
732-740.
[8] Durbeck, L., Macias, N., Weinsten, D., Johnson, C., andHollerback, J., SCIRun Haptic Display for Scientific
Visualization, Phantom Users Group Meetings, 1998.
[9] Enright, D., Marshner, S., and Fedkiw R., Animation and
Rendering of Complex Water Surfaces, ACM. Transaction
on Graphics (Proceedings of ACM SIGGRAPH2002), 21, 3,
2002, 737-744.
[10] Foster, N. and Fedkiw, R., Practical Animation of Liquids, In
Proceedings of SIGGRAPH2001, Annual Conference Series,
2001, 23-30.
[11] Hasegawa, S. and Sato, M., Real-time Rigid Body
Simulation for Haptic Interactions based on Contact Volume
of Polygonal Objects, Computer Graphics Forum(Proceedings of EUROGRAPHICS 2004), 23, 3, 2004, 529538.
[12] Hoerner, F. S., Fluid Dynamic Drag: Practical Information
on Aerodynamic Drag and Hydrodynamic Resistance,
Hoerner Fluid Dynamics, 1965.
[13] Kass, M. and Miller, G., Rapid, Stable Fluid Dynamics for
Computer Graphics, Computer Graphics (Proceedings of
ACM SIGGRAPH 1990), 24, 4, 1990, 49-57.
[14] Kim, S., Hasegawa, S., Koike, Y., and Sato, M., Tension
based 7-dof force feedback device: Spidar-g, In Proceedingsof IEEE Virtual Reality 2002, Annual Conference Series,
2002, 283-284.
[15]
Landau, L. D. and Lifschitz, E. M., Fluid Mechanics,
Pergamon Press, 1975.
[16] Mark, W. R., Randolph, S. C., Van Verth, J. M., and Taylor
II, R. M., Adding Force Feedback to Graphics Systems:
Issues and Solutions, In Proceedings of ACM SIGGRAPH 96,
Annual Conference Series, 1996, 447 - 452.
[17] Mcneely, W. A., Puterbaugh, K. D., and Troy, J. J., Six
Degree-of-Freedom Haptic Rendering Using Voxel
Sampling, In Proceedings of ACM SIGGRAPH 99, Annual
Conference Series, 1999, 401 - 408.
[18] Otaduy, M. A. and Lin, M. C., Sensation Preserving
Simplification for Haptic Rendering, ACM Transaction on
Graphics (Proceedings of ACM SIGGRAPH 2003), 23, 3,
2003, 543-553[19] Pai, D. K, Kees, van den D., James, D. L., Lang, J., Lloyd, J.
E., Richmond, J. L., Yau, S. H., Scanning PhysicalInteraction Behavior of 3D Objects, In Proceedings of
SIGGRAPH 2001, Annual Conference Series, 2001, 87-96.
[20] Ruspini, D. C., Kolarov, K., and Khatib, O., The Haptic
Display of Complex Graphical Environments, In
Proceedings of ACM SIGGRAPH 97, Annual Conference
Series, 1997, 345-352.
[21] Stam, J., Stable Fluids, In Proceedings of ACM SIGGRAPH
99, Annual Conference Series, 1999, 121-128.

[m/s]
0.4
0.2

vertical force

150
0

0
horizontal force-0.2

movement
direction’s force-0.4
0.5 1.0 1.5

2.0

2.5[s]

force

speed

[N]
300

speed

-150
-300

(a) t = 0.7(d) t = 1.3(c) t = 1.9(d) speed and force
Figure 8: Interaction of a sphere with water. Figs. (a)-(c) show different snapshot images during the interaction. The arrow
indicates the direction of the force due to water. (d) shows the relation between the speed of the sphere and the force. Threecomponents of the force are shown, that is, forces in the movement direction, the horizontal direction and vertical direction to
the movement
[m/s]
0.4
speed

0

-0.2

movement
direction’s force
-0.4
0.5

0

force

speed

0.2

[N]
horizontal 150
force
75

vertical -75
force
-150
1.0 [s](b) t = 0.75(d) speed and force(a) t = 0.6(c) t = 0.9
Figure 9: Interaction of a parallelopiped with water. The user moves the parallelopiped forward, however, the parallelopipedmoves left and right slightly due to the creation of wake vortices (see horizontal force in Fig. (d)).

PC

SPIDAR-G

stick(a) Setup for fish simulator.

(b) Virtual lure and rod.

(c) Example images during simulation.
Figure 10: A simulator for lure-fishing. (a) shows the haptic interface with a stick attached to simulate the rod. (b) shows thevirtual lure and the rod. (c) shows images during the simulation. In the upper left corner of the images in (c), the motion of theunderwater lure is displayed. The force on the lure and the force dragging the tip of the rod are indicated by the arrows.