Visual Simulation of Magnetic Fluids Using Dynamic Displacement Mapping for Spike
ShapesTomokazu ISHIKAWA†, Yonghao YUE†† (Member), Kei IWASAKI††† (Member),
Yoshinori DOBASHI†††† (Member), Tomoyuki NISHITA††††† (Member)†Tokyo University of Technology, †† Columbia University / JSPS, ††† Wakayama University / UEI Research†††† Hokkaido University / JST CREST, ††††† UEI Research / Hiroshima Shudo University<Summary> Although several simulation methods have been proposed to analyze the behavior of magnetic fluids in the computational physics field, it is still not reasonably tractable to simulate the spiking phenomenon using fully physically based
methods. To synthesize spike shapes, a procedural approach has been proposed in the computer graphics field recently for visual simulation, which, however, cannot simulate dynamic arrangement and deformation of spikes since the spike shapes are determined beforehand. To overcome this drawback, we propose an improved model to incorporate the dynamic behavior of the
spikes. The model is based on the assumption that the spikes are attracted by an external magnetic field and that each spike repels the others because of their magnetization. The position of the spikes on the liquid surface are determined in the simulation,and the spike shape is mapped according to the direction of the magnetic field lines. The coordinates at which the spike protrudes from the liquid surface are calculated in the fluid simulation.
Keywords: magnetic fluid, fluid simulation

1.

Introduction

In the field of computer graphics, the visual simulation of fluid

Solvent

phenomena is an important research topic. Several methods have

(Water or Oil)been proposed for representing realistic incompressible fluids bysolving the Navier–Stokes equations. In some fluid phenomena, theSurfactant
Magnetic microparticle

interaction between a magnetic field and a fluid is important: forexample the aurora borealis 1) is an interplay between the fluid
behavior of plasma and the magnetic field of the earth. In this paper,
we focus on the interaction between magnetic fluids and theFig. 1 Composition of magnetic fluid

external magnetic field, and propose a method for representing theresulting fluid behavior. The difference between a magnetic fluid
and the aurora is the presence and absence of a free boundary
surface, respectively.
A magnetic fluid is a colloidal solution consisting of
microparticles of ferromagnetic bodies, a surfactant that covers theFig. 2 Spiking phenomenon of magnetic fluids (photograph)magnetic microparticles, and a solvent that acts as the matrix (seeFig. 1). When a magnet is located near a magnetic fluid, the surface

proposed a method for synthesizing these spike shapes byof the fluid forms spikes (such as horns) along the direction of theemploying a procedural approach 2). Since the spike shapes weremagnetic field generated by the magnet. This is known as thedetermined beforehand, this method could not simulate the“spiking phenomenon” (see Fig. 2). This phenomenon is used to

dynamic arrangement and deformation of the spikes. In this paper,

create artworks, where interesting shapes are generated by applying

we propose a model to incorporate such dynamic arrangement andmagnetic forces to the fluids.

deformation of the spikes.

Within the field of computer graphics (CG), Ishikawa et al.

1

Magnetic field line
Spike

Attractive force tothe center

Seed
Liquid surface

Particle
External
net

Spike

Repulsive force
between the spikes

magMagnetic field line

External
magnetFig. 4 Forces acting on “seed”

Fig. 3 Spike and magnetization of particles

SPH particles). Each particle represents a small magnetic fluid

2. Related Workelement including microparticles, surfactant and a solvent, and its
motion is calculated by taking into account the properties of both

Many methods have been proposed to simulate incompressible

fluid and magnetic body. After calculating the fluid surface using

fluids such as smoke and flames 3), 4). Goktekin et al. proposed a

SPH particles, we map each spike shape to the liquid surface. Un-

simulation method for viscoelastic fluids by incorporating an elastic

like the previous method 2), we model each spike individually,5)

term into the Navier–Stokes equations . Stam and Fuime intro-

hence allowing the spikes to rearrange dynamically. When using

duced the smoothed-particle hydrodynamics (SPH) method into thethe SPH method, we place “markers” on the fluid surface for map-

CG field for representing flames and smoke 6). Müller et al. pro-ping the spike shapes. Each spike shape is represented as a super-

posed an SPH-based method to simulate fluids with free surfaces 7).

position of trigonometric functions. We need to calculate only the

A method for calculating the density distribution as voxel data iscenter of the mapping (i.e., the place of the marker) which we de-typically used for the fluid modeling of phenomena such as thenote as a “seed.”aurora borealis. However, fluid phenomena such as water which

By calculating the movement of the seeds, considering the fluid

have free boundary surface are often simulated using the particle

behavior and the force between seeds, we map the spike shapes tomethod. Our method also uses an SPH-based method 2) to simulate

the seeds at each time step of the simulation. The SPH particles andmagnetic fluids.

seeds are simultaneously computed. The number of seeds consid-

In the field of physics, the characteristics of magnetic fluids haveered for the calculations is proportional to the external magnetic

been studied since 1960. Han et al. modeled the formation of aforce.

chain shape between colloidal particles according to the magnetization of the particles 8). Combined with the lattice Boltzmann meth-

3.1. Movement of “seeds”

od, they showed that colloidal particles would form lines along the

As a seed is a mapping point, it should be located on the liquidmagnetic field. However, their method cannot accurately reproduce

surface. In our method, a seed is modeled as a special type of SPH

spike shapes. Yoshikawa et al. combined the moving-particle

particle, whose mass is lower than that of a standard SPH particle.

semi-implicit (MPS) method with the finite element method (FEM)

That is, mass of seed particles is m if mass of standard SPH parti-

and simulated magnetic fluids 9). Even when using 100,000 parti-

cle is m. Let  be 0.9 in this paper. At a seed location, the pressure

cles and a mesh with 250,000 tetrahedra, they were able to repro-

from other SPH particles is calculated. However, the pressure fromduce only a single spike.

the seeds to the other SPH particles is not considered, since these

We observed that the number of spikes increases when a strong

seeds are not physically based entities for representing the fluid.

magnetic force is applied to magnetic fluids, and the spikes extend

Note that the seeds are considered in the density calculation. A bar

along magnetic field lines. It seems these dynamic behaviors aremagnet is placed outside the fluid; we call this magnet an external

not reasonably tractable using fully physically based method.

magnet, which is used to magnetize the fluid particles. The fluid

Hence, we instead propose a physically plausible method to repro-

particles are hence magnetized due to this external magnet. In eachduce these dynamic features.

spike, the (real) magnetized micro particles are arranged along themagnetic field line and oriented by attractive forces (Fig. 3).

3. Our Model

Therefore, we assume that the spike itself forms a single magnet,
thus explaining why separate spikes repel each other. In addition,

Our method is based on the SPH method, which is used to calcu-

attractive forces act toward the intersection point between thelate the magnetic forces acting on simulation particles (hereinafter2

Table 1: Parameter setting used in equationsuv

param.

meaningvalue

dt

timestep

0.01 sporder of the attractive force term

2.0qorder of the repulsive force term

4.0



depth of the potential well

0.02 J

collision diameter2.6×10−3 m

proportional coefficient in Eq. (4)15.0×10−3 m

Refer to only this spikeymethod, we can determine the coordinates of the seed at which the
spike is formed.

x3.2. Mapping the spike shapes
Mapping
After calculating the coordinates of the seeds, we map the spikearea

shape. The method is based on displacement mapping, where the
Displacement along

fluid surface of the seed periphery is directly deformed. In addition,

the magnetic field line
z

Spike
w

we allow the fluid surface to displace in the direction of the magnetic field line (because the spike of the magnetic fluids protrudes

Liquid surfacexin this direction), as indicated by the height map (Fig. 5 bottom).
The surface shape for mapping is taken from an analysis of the casein which a vertical magnetic field is applied to the computation of

Seed
Fig. 5 Spike shape reproduced by a function (top), Conceptual
diagram of mapping to the liquid surface (middle: tophorizontal surface of a magnetic fluid 10). The function is known as
the one that we can analytically seek for. This shape can be represented as a height field as follows:

view, bottom: cross-sectional view)kw  u, v   C0  cos2 3u  v   cos k2  3u  v   cos kv  .

(3)

liquid surface and the center of the magnetic field line induced bywhere u, v, and w are the three dimensions. C0 indicates the ampli-

an external magnet, as the top of the external magnet is atude of spike and k represents the frequency. We assume that C0 is

well-balanced point (Fig. 4).

proportional to the magnitude of the magnetic field as follows,

We propose that the stable point of the spike can be determined inrelation to the following two forces: the attractive force toward theC0   H  r  ,(4)center and the repulsive force between the spikes. In this method,where  is the proportional coefficient, and |H(r)| is the magnitude

we calculate the mapping points (seeds) by calculating the interac-of the magnetic field at seed position r. H(r) is determined by thetion of the forces between the seeds, rather than those between theposition of an external magnet and the magnetization of each SPH

spikes. We consider the Lennard-Jones potential (Eq. (1)), where dparticle. Please refer to the Appendix for the calculation of thisis the distance between any two seeds, and p, q, , and  are pa-

magnetic vector. We map each of the spikes by calculating Equa-

rameters. The attractive force and repulsive force are included intion (3), referring to the height map within a certain radius of thethe Lennard-Jones potential.

center of the top of the spike (Fig. 5 top). A radius that contains a   p   q U  d   4       .
dd    single spike has length equal to half of the cycle length in Equation(1)(3). We create a spike shape by selecting the shortest distance between the seeds as the wavelength. When the shortest distance is

The Lennard-Jones potential is typically used for computing thedmin, the parameter k of the frequency in Eq. (3) is determined asinteraction between two atoms. The force F acting on the seed onfollows,

account of this potential is calculated by the following formula.
F  d   U  d  .

k(2)

By calculating the force from the distance between the seeds and4 3
.
3d min(5)

At each seed point, we map the spike shape to the seed position
corresponding to the origin of the height map (Fig. 5 middle).

considering their movement due to an external force in the SPH
3

SPH particles representingzmagnetic fluid
y

x

Seed
(a) Initial state(b) Magnet placement for Figure 7(c) Magnets placement for Figure 9

Fig. 6 The simulation space for the proposed method, (a) The magnetic fluids, which are represented by a set of SPH
particles, are stored in a cubic container, (b) The magnet is placed underneath the container, (c) The magnets areplaced in a setup so as to recreate the artwork in Fig. 9(a) 300 timestep(b) 330 timestep(c) 360 timestep

Fig. 7 Resulting images from the proposed method, showing simulation results when the magnet approaches from thebottom of the magnetic fluid

Fig. 8 Comparison of real spike phenomena (top), previous method (middle) and the proposed method (bottom)

4. ResultsCore 2 Duo™ 3.33GHz, RAM: 3.25GB, GPU: NVIDIA GeForceGTX 480). The surfaces of the magnetic fluids were extracted us-

For the simulation, we used C++ and OpenGL for CPU program-ing the method proposed by Yu et al. 11). Using the above men-ing, a CUDA-based GPU for the SPH method and the calculation

tioned PC, the average computation time for the surface construc-of the motion of the seeds and SPH particles, and POV-Ray fortion of each frame was 2 min. The parameters used in our simula-rendering. The average computation time of the simulation for ation are shown in Table 1.

single time step was 0.16 s when using a standard PC (CPU: Intel®
4(a) 350 timestep(b) 390 timestep

Fig. 9 Simulation results emulating an artwork, we set one magnet above and four magnets below the magnetic fluid
The initial fluid surface is shown in Fig. 6 (a). The magnetic fluids,
which are represented by a set of SPH particles including the seeds,
are stored in a cubic container. To achieve simulation results suchas those shown in Figs. 7 and 8, the magnet is placed underneaththe container. To emulate the art work in Fig. 9, we set the magnets
as shown in Fig. 6 (c).
(a) 370 timestep

Fig. 8 shows the comparisons between real magnetic fluids, theprevious method and proposed method by changing the magnetic
force. The images in the middle row of images in Fig. 8 show theresult using the previous method 2). The previous method has a
restriction in which the number of spikes is fixed, hence cannot
represent varying number of spikes. In addition, the arrangement of(b) 390 timestep

the spikes cannot be dynamically changed in this method. We canreproduce that the number of spikes increases when a strong magnetic force is applied to the magnetic fluids, as well as that the dynamic arrangement of spikes.
Fig. 9 shows the results to emulate an art work 12), 13).
Fig. 10 shows simulation result when the magnet goes away fromthe bottom of the magnetic fluid. The proposed method can achieve(c) 410 timestep

a wide range of expression since our method can change the loca-(d) 430 timestep(e) 450 timestep
Fig. 11 The strokes show the magnetic field lines, we
Fig 10 Simulated images showing the case where themagnet is moved away from the bottom of thecan see that the spikes extend along the magnetic field lines

magnetic fluid and the spikes disappear57) M. M¨uller, D. Charypar, M. Gross: “Particle-Based Fluid Simulation fortion of spikes dynamically.
Fig. 11 shows the result images of spike shape with magnetic field

Interactive Applications”, Proc. of the 2003 ACM SIGGRAPH/ Eu-lines. We can see that each spike extends along the magnetic fieldrographics Symposium on Computer Animation, pp. 154–159 (2003).
8) K. Han, Y. T. Feng, D. R. J. Owen: “Three-Dimensional Modeling andlines, like actual magnetic fluids.

5. Conclusions and Future Work
By combining the SPH method with moving the mapping point of

Simulation of Magnetorheological Fluids”, International Journal for
Numerical Methods in Engineering, 2010, Vol.84, No.11, pp. 1273–1302(2010).

spikes, we have represented the dynamic deformation of spikes in9) G. Yoshikawa, K. Hirata, F. Miyasaka, Y. Okaue: “Numerical Analysis ofmagnetic fluids. These spikes were successfully incorporated byTransitional Behavior of Ferrofluid Employing MPS Method and FEM”,

incorporating the calculation of the magnetic field into the SPH

IEEE Transactions on Magnetics, Vol.47, No.5, pp. 1370–1373 (2011).

method. By moving mapping points along the fluid flow, we calcu-

10) R. Rosensweig: “Magnetic Fluids”, Annual Review of Fluid Mechanics,

lated the coordinates for mapping the spike shape.

Vol.19, pp. 437–461 (1987).

In the future study, we would like to include other factors that de-11) J. Yu, G. Turk: “Reconstructing Surfaces of Particle-Based Fluids Using

pend on the magnetic field, such as surface tension and magnetic

Anisotropic Kernels”, Proc. of the 2010 ACM SIGGRAPH/ Eu-

susceptibility. Moreover, we are planning to make interactive ap-

rographics Symposium on Computer Animation, pp. 217–225 (2010).

plications for studying magnetic fluids.

References

12) S. Kodama: Protrude, Flow, In Art Gallery of SIGGRAPH 2001 (2001).
13) http://sachikokodama.com/ .

Appendix

1) T. Ishikawa, Y. Yue, K. Iwasaki, Y. Dobashi, T. Nishita: “Modeling of
Aurora Borealis Using the Observed Data”, Proc. of SCCG 2011, pp.

A. SPH method considering the magnetic force

35–38 (2011).
2) T. Ishikawa, Y. Yue, K. Iwasaki, Y. Dobashi, T. Nishita: “Visual Simula-

The behavior of incompressible fluids is described by the follow-

tion of Magnetic Fluids”, Proc. of GRAPP 2012, pp. 319–327 (2012).

ing equations.

3) J. Stam, “Stable Fluids”, Proc. of SIGGRAPH 1999, pp. 121–128

  u  0,

(A1)

u
1
   u   u  p  2u  F.
t


(A2)

(1999).
4) R. Fedkiw, J. Stam, H. W. Jensen: “Visual Simulation of Smoke”, Proc.
of SIGGRAPH 2001, pp. 15–22 (2001).

Eq. (A1) is the continuity equation and Eq. (A2) (Navier-Stokes

5) T. G. Goktekin, A. W. Bargteil, J. F. O ’Brien: “A method for Animating

equation) describes the conservation of momentum. u is the veloc-

Viscoelastic Fluids”, Proc. of SIGGRAPH 2004, pp. 463–468 (2004).

ity vector, t is time,  is the fluid density, p is the pressure and  is6) J. Stam, E. Fiume: “Depicting Fire and Other Gaseous Phenomena

the kinematic viscosity coefficient. F is the external force that in-Using Diffusion Processes”, Proc. of SIGGRAPH 1995, pp. 129–136

cludes the gravity, the magnetic force, and the surface tension. We(1995).

solve these equations by using SPH method.

(a)

(b)

Fig. A Calculation of the magnetization and the magnetic force, (a) First, we calculate the magnetic field vector at each SPH particle
induced by a magnetic dipole, (b) Next, we calculate the influence of other particles from the magnetized particles
6Each SPH particle represents a small magnetic fluid element, andTomokazu Ishikawa

motions of the SPH particle and the mapping point “seed” are cal-Tomokazu Ishikawa is a Research Associate in
School of Media Science at Tokyo University ofTechnology. He received the B.S. degrees inInformation Science from Tokyo University of
Science in 2003, and his M.S. and Ph.D inComplexity Science and Engineering in 2005,
2012, respectively. He was hired as a systemengineer by IBM Japan, Ltd. until 2007. His
research interests center on computer graphics,
including visual simulation.

culated by taking into account the magnetic force. By applying an
external magnet, SPH particles are magnetized and work as small
spherical magnet. Therefore, our method first computes the magnetic field at equilibrium state, taking into account the magnetization of SPH particles. After that, the magnetic force for each SPH
particle is calculated. We call the magnetic vector induced by the

Yonghao Yue

netic vector field. Thus, in order to obtain the final magnetic vector(Member)
Yonghao Yue is currently a postdoctoral research
fellow at Columbia University with the fellowship support by Japan Society for the Promotionof Science (Postdoctoral Fellowship for Research Abroad). He received his Ph.D. degree inInformation Science and Technology from the
University of Tokyo in 2011. He was an assistantprofessor in the Department of Complexity
Science and Engineering at University of Tokyo
from 2011 to 2013, and an assistant professor in
the Department of Information Science at University of Tokyo from 2012 to 2013. His research interests include physical based renderingand simulation.

H(rj) at particle j, the magnetic interactions between particles have

Kei Iwasaki

magnetic bar as a background magnetic vector. Let us assume thatthe origin is at the midpoint between the north and the south poles.
Then, the background magnetic vector Hdipole(r) at position r isexpressed by the following equation.
H dipole  r   

1
m r
 3 ,
4
r

(A3)

where μ is the permeability of the magnetic fluid, m is the magnetic
moment and r = |r|. Each SPH particle is magnetized due to thebackground magnetic vector field and induces an additional mag-to be computed by solving the following equation,H  r j   H dipole  r j  

V N  H  r j   rij
 r 3 ,4 i 1
ij

(A4)

i j

(Member)
Kei Iwasaki is currently an associate professor in
the Faculty of Systems Engineering at Wakayama University. He received his Ph.D degreefrom The University of Tokyo in 2004. His
research interests include real-time renderingand visual simulation.

where V is the volume of a particle and we assume that the volumeof all particles are equal, ri is the position of particle i, N is the total
number of particles,  is the magnetic susceptibility. rij = rj − ri, andrij = | rj − ri |. We calculate Eq. (A4) by taking into account the
influences from all the particles. We use the Gauss-Seidel method
to solve Eq. (A4).
Finally, we obtain the magnetic force based on the distribution of
the magnetic field by the following equation 10),
Fmag  ri   

 H  ri 
2

2

.

(A5)

（Received July 2, 2013）（Revised Nov. 2, 2013）

Yoshinori Dobashi(Member)
Yoshinori Dobashi is an associate professor at
Hokkaido University in the Graduate School ofInformation Science and Technology, Japan
since 2006. His research interests center oncomputer graphics including realistic image
synthesis, efficient rendering, and simulation ofnatural phenomena. Dobashi received his BE,
ME and Ph.D in Engineering in 1992, 1994, and1997, respectively, from Hiroshima University.
He worked at Hiroshima City University from1997 to 2000 as a research associate.

Tomoyuki Nishita(Member)
Tomoyuki Nishita is a professor in the Department of Economical Sience, Hiroshima Shudo
University, and was in the Department of Complexity Science and Engineering (also in the
Department of Information Science) at University of Tokyo, Japan since 1998. He taught at
Fukuyama University from 1979 to 1998. He
was an associate researcher in the Engineering
Computer Graphics Laboratory at Brigham
Young University from 1988 to 1989. His research interests center in computer graphics
including lighting models, hidden-surface removal, and antialiasing. Nishita received his BE,
ME and Ph.D in Engineering in 1971, 1973, and1985, respectively, from Hiroshima University.
Dr. Nishita received Research Award on Computer Graphics from Information Processing
Society of Japan in 1987, and also received
Steaven A. Coons award from SIGGRAPH in
2005. He was a president of The Institute of
Image Electronics Engineers of Japan in 2009.

7

8