A Current Loop Model for the Fast Simulation of Ferrofluids

Ferrofluids are oil-based liquids containing magnetic particles that interact with magnetic fields without solidifying. Leveraging the exploration of new applications of these promising materials (such as in optics, medicine and engineering) requires high fidelity modeling and simulation capabilities in order to accurately explore ferrofluids in silico. While recent work addressed the macroscopic simulation of large-scale ferrofluids using smoothed-particle hydrodynamics (SPH), such simulations are computationally expensive. In their work, the Kelvin force model has been used to calculate interactions between different SPH particles. The application of this model results in a force pointing outwards with respect to the fluid surface causing significant levitation problems. This drawback limits the application of more advanced and efficient SPH frameworks such as divergence-free SPH (DFSPH) or implicit incompressible SPH (IISPH). In this contribution, we propose a current loop magnetic force model which enables the fast macroscopic simulation of ferrofluids. Our new force model results in a force term pointing inwards allowing for more stable and fast simulations of ferrofluids using DFSPH and IISPH.


INTRODUCTION
T HE simulation of complex natural phenomena is an active topic of research in several scientific communities ranging from computational physics to visual computing.At least since Terzopoulos et al.'s pioneering work [1] on elastically deformable models presented at SIGGRAPH 1987, physically-based simulation is also at the core of the computer graphics community.Twelve years later, Stam [2] presented his seminal work on stable fluids establishing fluid dynamics research in computer graphics.Since then graphics researchers have managed to continuously push the boundaries in fluid simulation such as realistically simulating water [3] and smoke [4] as well as fluid flow around objects [5].Recent contributions, e.g., address fluid phenomena such as bubble rings [6] and waves on a vast ocean [7], and even complex natural phenomena such as storms [8], several weather effects [9], and wildfires [10].In 2007, Bridson and M€ uller-Fischer [11] gave an influential course at SIGGRAPH providing a practical introduction to fluid simulation for graphics enabling the audience to animate fully three-dimensional incompressible flow.Significant research effort has been devoted to enable accurate and efficient fluid simulation.Among others, the powerful particle-based (Lagrangian) smoothed particle hydrodynamics (SPH) technique has been introduced to the graphics community for simulating fluids by M€ uller-Fischer et al. [12].Previously, Desbrun and Gascuel [13] already simulated highly deformable bodies using SPH.Koschier et al. [14] later published an excellent Eurographics tutorial focusing on state-of-the-art fluid simulations based on SPH.Moreover, several grid-based (Eulerian) approaches [15] became popular as well as hybrid Lagrangian-Eulerian approaches such as the fluid-implicit-particle (FLIP) method [16].Aiming for more realism, researchers have incorporated several physical models, e.g., for surface tension [17], viscosity [18], [19], and turbulent effects [20] into their fluid solvers.
Inspired by the interesting behavior and practical relevance of so-called ferrofluids (i.e.oil-based liquids which contain magnetic particles), the graphics community has recently addressed their simulation focusing at the macroscopic scale.In the absence of an external magnetic field, these materials behave like regular Newtonian fluids.Once a magnetic field is applied, interesting complex patterns are generated showing characteristic spikes oriented according to the field lines of the external magnetic field.Covering the underlying physics and -at the same time -handling the complicated underlying surface geometry fits the skill-set of the computer graphics community quite well.Huang et al. [21] from this community have been the first to address the first-principle-based macroscopic simulation of ferrofluids.Huang et al. utilized the SPH methodology and incorporated the appropriate handling of the electromagnetic effects.Later, Huang and Michels [22] devised a boundary element method for the simulation of ferrofluids and Ni et al. [23] introduced a level-set method for magnetic Substance Simulation.These contributions aim for facilitating the exploration of new applications of these promising materials such as in optics, medicine and engineering by providing modeling and simulation capabilities in order to accurately explore ferrofluids in silico.In the conmathrm of electromagnetism, graphics researchers have previously worked, e.g., on magnetic rigid bodies addressing magnets in motion [24] and magnetization dynamics for magnetic object interactions [25].While Huang et al. have proposed an approach for the macroscopic simulation of ferrofluids, their Kelvin magnetic force model leads to an outward-pointing magnetic force, which would levitate the particles around the fluid surface and leads to an unphysical surface shape.In terms of SPH framework, Huang et al. have utilized the weakly compressible SPH (WCSPH) in contrast to state-of-the-art SPH methodology such as divergence-free SPH (DFSPH [26]) or implicit incompressible SPH (IISPH [27]).Compared to DFSPH and IISPH that only penalize fluid compression, WCSPH also penalizes fluid expansion which happens on the surface of previous ferrofluid simulations using the Kelvin force model.This could alleviate the levitation problem mentioned above.However, significant instabilities are observed when using WCSPH with large temporal time steps resulting in low performance.In this contribution we aim to overcome this limitation by introducing a new magnetic force model enabling the fast simulation of ferrofluids.Our so-called current loop magnetic force model results in a force term pointing inwards on the fluid surface, which exempts the requirement to penalize fluid expansion, and thus can be well integrated with DFSPH and IISPH.A large temporal step size can be used in our simulations without instability issues.Moreover, we utilize a better linear solver for the computation of the magnetic field.We demonstrate the capabilities of our new approach by presenting various numerical experiments and a comparison to Huang et al.'s previous model.

RELATED WORK
In the following, we discuss related work on fluid dynamics followed by previous work on ferrofluids published within the computer graphics community.We would like to point out that the simulation of ferrofluids is also an active topic of research in different scientific communities [28].However, macroscopic simulations of ferrofluids covering a variety of visually recognizable aspects can almost exclusively be found within computer graphics research.

Fluid Simulation
In general, there are two popular types of approaches in fluid dynamics based on Eulerian and Lagrangian representations.Following the Eulerian way, the fluid flow is represented by a velocity field.The flow velocity is given by uðx; tÞ defined at a certain position x on a grid at time t.In contrast, following the Lagrangian way, particles are introduced representing individual fluid elements with a flow velocity uði; tÞ attached to each particle i.It is relatively easy to discretize the underlying differential equations using the Eulerian approach.However, dealing with meshes -especially with larger shape variations or topology changesintroduces a significant complexity overhead.Within the graphics community, Lagrangian approaches such as SPH [12] and hybrid approaches such as FLIP [16] are popular choices.
SPH utilizes kernels to estimate the physical quantities at the particle positions.The corresponding governing equations are discretized with respect to each particle.The first work on SPH dates back to 1977 [29], [30] focusing on problems in astronomy before Desbrun and Gascuel [13], and M€ uller-Fischer et al. [12] made SPH popular in graphics.
To deal with incompressibility, weakly compressible SPH (WCSPH [31], [32]) utilizes stiff equations of state (EOS) regarding pressure and density.To ensure stability and avoid oscillations, the CFL condition must be satisfied and thus a small time step is required in conjunction with explicit time integration, which would greatly increase the computational cost.Solenthaler et al. proposed predictive-corrective incompressible SPH (PCISPH [33]) which achieves incompressibility while allowing for larger time steps.With superior scaling properties, implicit incompressible SPH (IISPH [27]) achieves speedups over PCISPH.Divergencefree SPH (DFSPH [26]) further enforces the incompressibility on both position and velocity levels.This ensures realistic behavior and results in a stability increase.For boundary handling, Akinci et al. [34] proposed a two-way coupling approach.Adami et al. [35] developed a wall model which prevents particle penetrations based on a local force balance.Koschier et al. [36] used a density map to avoid the sampling of boundary particles.Bender et al. [37] further developed volume maps to reduce computation time and memory requirements.Regarding surface tension, ghost SPH [38] has been proposed to deal with surface tension dispersions.Akinci et al. [39] proposed a surface tension and adhesion force model which can handle large surface tension.He et al. [40] proposed a surface tension model based on a free surface energy functional.Zorilla et al. [17] accelerated the computation of surface tension by classifying the particle and Monte Carlo integration.
The early work of hybrid Lagrangian-Eulerian approaches dates back to the particle-in-cell (PIC) approach [41].To address the problem of large numerical dissipation using PIC, the fluid-implicit-particle (FLIP) method was developed by Brackbill et al. [42].Within the graphics community, FLIP was introduced by Zhu et al. [16].Batty et al. [43] proposed an efficient and robust method to handle the irregular boundary geometry.Ng et al. [44] further improved the fluid-solid coupling at the boundaries which converges in the L 1 -norm.Boyd et al. [45] developed a multi-phase FLIP method for accurate surface tension simulation.Ferstl et al. [46] introduced narrow band FLIP to reduce particle counts and computation time.

Ferrofluids
Early research results on ferrofluids can be found within the physics community [47], [48].Several contributions focused on their simulation [49], [50], [51], [52] covering properties such as the formation of spikes.However, these simulations cannot handle complex dynamics.Instead, the process towards the equilibrium state is simulated.Liu et al. [53] simulated ferrofluid droplets and Zhu et al. [54] utilized a levelset method for the simulation of ferrofluids.
Within the graphics community, Ishikawa et al. [55], [56] used SPH to simulate ferrofluids as particles.However, the spiky fluid surface is generated by procedural modeling.Huang et al. [21] have been the first who devised an SPH approach to the qualitatively accurate simulation of large-scale ferrofluids based on physical principles.In contrast to previous work, the characteristic spike patterns of ferrofluids are generated in a physical meaningful way without a priori knowledge of the spike pattern.The Kelvin force model is employed and incorporated into the SPH fluid solver.Later, Ni et al. [23] utilized a level-set method to simulate magnetic substances.In their work, the magnetic force is modeled as the interfacial Helmholtz force drawn from the Minkowski form of the Maxwell stress tensor.Huang and Michels [22] devised a surface-only approach which uses a boundary element method to solve the magnetization problem and compute the magnetic pressure for force evaluation.Sun et al. [57] utilized the material point method to simulate nonlinearly magnetized materials.
On a different trajectory, the simulation of magnetic rigid bodies has also been addressed in the graphics community.In this regard, Thomaszewski et al. [24] simulated the rigid magnet by calculating the interaction between the subdivided magnetic dipole and the external field.Kim et al. [25] improved the simulation based on magnetization dynamics which enabled the incorporation of mutual inductance and remanent magnetization.Kim and Han [58] further increased the stability and enabled the simulation of magnets with arbitrary shape.

METHODOLOGY
In this section, we will introduce our novel approach.First, the general information of our fluid solver will be discussed.This includes the control equations, the SPH framework DFSPH [26], the viscosity force and the surface tension force models, and the boundary handling.After that, the computation of the magnetic field is discussed which is mainly a recap of previous work from Huang et al. [21].This is followed by the introduction of our novel current loop force model.We derive its new formulation and illustrate how this leads to more stable simulations of ferrofluids.

Fluid Solver
The temporal evolution of the fluid's state can be simulated by solving the Navier-Stokes equations.Incompressible fluid is described by the continuity equation, and the momentum equation, in which the density is denoted by r, the pressure by p, the viscosity coefficient by m, and f ext represents the external net force which usually contains the magnetic forces as well as gravity and surface tension.Please note, that the vector Laplacian is similarly defined as its scalar counterpart and simply acts component-wise.SPH model discretizes the fluid domain using a set of particles.A physical quantity A at position r i is approximated by summation of its neighbouring values in which W ij ¼ W ðr i À r j Þ and W ðrÞ refer to a kernel function [30].
The density is estimated by which is a weighted summation of the particle masses.Since the mass is equally distributed among the particles and a constant value, the continuity Equation ( 1) is naturally satisfied.The left hand side of the momentum Equation ( 2) can be re-written into a material derivative form r Dv=Dt.In a SPH framework, pressure and viscosity are modeled as forces acting on the particles.Following a certain particle i, we have the control equation with the force while g denotes the gravity.In a standard SPH [12] or WCSPH framework [31], the pressure force is computed using the equation of state.However, to ensure incompressibility, high stiffness is required which limits the step size used for the numerical integration.In IISPH [27] and DFSPH [26], iterative methods are proposed to ensure incompressibility while maintaining a large time step size.In IISPH, the pressure is computed iteratively by solving a linear system.In DFSPH, a constant density solver is used to eliminate the density error and a divergence-free solver is used to ensure the divergence-free velocity field.
In our work, we utilize both IISPH and DFSPH for comparison experiments and mainly use the DFSPH framework to conduct numerical experiments.Algorithm 1 gives an overview of DFSPH.Please note, that a is a factor which is depending on the particle positions.For details about the constant density solver and the divergence-free solver, please refer to Bender et al. [26].
The standard viscosity force is modeled as f viscosity i ¼ nr 2 v i which requires the SPH discretization of the Laplacian operator However, the use of second derivatives is problematic because they cannot capture the differences of physical quantities and lead to noise.We employ the viscosity model from Weiler et al. [18] which uses a combination of first derivative and finite differences with dimension number d. Please note, that the 10 À2 h 2 term in the denominator is used to avoid the singularity.We follow the work of Akinci et al. [39] to model surface tension.In their work, both the cohesion effect between particles and the surface area minimization are addressed The cohesion term is given by in which s denotes the surface tension coefficient and C denotes a spline function.The surface minimization term is given by in which n i provides normal information avoiding the explicit computation of the curvature which would require second-order derivatives.n i and n j are not normalized such that they are close to zero inside the fluid domain, which leads to a surface tension force close to zero.We refer the readers to Akinci et al. [39] for more discussions about the surface tension model.For boundary handling, we follow the work of Bender et al. [37] which precomputes a volume map.The volume map provides an implicit boundary representation.It utilizes the narrow band data structure to store the information and query values at run-time.

Magnetic Field
Following the previous section, the remaining work towards the simulation of ferrofluids is a proper magnetic force model.The first step is to calculate the magnetization when an external magnetic field is imposed on the fluid particles.This part mainly follows the work from Huang et al. [21].
The magnetic field is usually described by the magnetic flux density B and the magnetic field strength H.The relationship between these quantities is given by where the magnetization vector field is denoted by M and vacuum permeability is denoted by m 0 .Please note, that MðrÞ is a field vector which describes the magnetic moment's density at location r.In our ferrofluid simulation, each fluid particle is naively equated to a tiny magnet, thus the magnetic moment m is used to describe the magnetic strength of a single particle.In the SPH framework, the relationship between the magnetization vector field M and the particle magnetic moment m i is given by where the kernel function is denoted by W ðr À r i ; hÞ.Under a certain external magnetic field H ext , the ferrofluid will be magnetized with magnetic moment m i for each particle to be determined.These magnetized particles would generate a magnetic field H ferro .Thus According to Maxwell's equation, the magnetic flux density B is divergence-free and r Á H ext ¼ 0 outside the external magnet .Following the work of Huang et al. [21], H ferro is curl-free and it can be represented with H ferro ¼ Àrf, where f is a scalar function.Putting this back into Eq.( 15), we end up with a Poisson equation Omitting the derivation, we directly use the solution to the Poisson's equation above from Huang et al. [21] This equation describes the magnetic field at position r generated by a single ferrofluid particle centered at r ¼ 0 with magnetic moment m i , where r ¼ jrj, r ¼ r=r, and W avr ðrÞ is the averaged kernel value contained within a ball of radius r The magnetic field generated by the ferrofluid is a summation of contributions from individual particles Substituting Eqs. ( 13) and (20) into Eq.( 14) regarding each particle position r i , we obtain The above equation contains two unknowns Bðr i Þ and m i .To deal with this, we will use the approximation with volume V ¼ dx 3 .Please note, that the particles are initially positioned on a uniform grid with spacing dx.Moreover, the constitutive relation is given by where x is the magnetic susceptibility, which describes the level of magnetization of a substance in response to an imposed magnetic field.Combining Eqs. ( 12), ( 22) and ( 23), we obtain a linear relationship between Bðr i Þ and m i Substituting this back into Eq.( 21), we obtain a linear system where b i ¼ Bðr i Þ=m 0 corresponds to the magnetic flux density at particle position r i , and G ¼ H ferro þ M is the summation operator in Eq. ( 21) regarding the particle magnetic moment m i .At each time step, the particle position r i and the external field H ext ðr i Þ are given.We solve Eq. ( 25) for b i and use Eq. ( 24) to get the particle magnetic moment m i .
In the previous work of Huang et al. [21], Eq. ( 25) is solved using least squares conjugate gradient.However, we found that a standard conjugate gradient method is sufficient.

Current Loop Force Model
After the magnetization process, the magnetic moments are usually aligned with the magnetic field, so that we can ignore the torque and only consider the force imposed on the ferrofluid particle.
According to the theoretical work of Byrne [47], there are multiple magnetic force models, and in the work of Huang et al. [21] the Kelvin force model is adopted with where ðA Á rÞB ¼ ð@B i =@x j ÞA j .
In this work, we use another magnetic force model: a current loop model with the following formulation: where ðrBÞ Á A ¼ A i ð@B i =@x j Þ.If both, A and B, are given as column vectors, then ðrBÞ Á A ¼ ðrBÞ T A which returns a column vector resulting from multiplying a matrix with a column vector.According to Byrne [47], the first term m 0 ðrHÞ Á M in Eq. ( 27) is equivalent to the Kelvin force model of Eq. (26).Define F M ¼ m 0 ðrMÞ Á M, then Eq. ( 27) can be rewritten as Hence, our current loop force model is the result of adding a bulk term to the Kelvin force model.The Kelvin model and the current loop model have different force densities, but when combined with corresponding surface traction discontinuity, they have the same total pressure jump on the surface induced by magnetic forces as explained in the conclusion of Byrne [47]: "When [many classical theories including Kelvin and current loop forces are] augmented by the corresponding differential surface traction to account for boundary layer forces, a unique formulation is obtained for the pressure rise to a containing boundary, and a unique surface integral method for total body force."In the continuous form, the equilibrium shape of a ferrofluid surface is determined by the pressure discontinuity across the ferrofluid surface.This is analogous to the different behavior of water and mercury in capillary effects.The concave and convex shape of the liquid surface in the thin tube has opposite capillary action, which is decided by the pressure jump.The current loop model gives an inward-pointing magnetic force which leads to more stable simulation as we will show later.
In the discretized form, we first consider the magnetic force from a source particle towards a target particle The first integral is given by Huang et al. [ The second part on the right side of Eq. ( 29) is addressed as follows: where and in the above W 0 ðr; hÞ is the derivative of the kernel function W ðr; hÞ.Hence, we obtain Z F M dr ¼ m 0 Z W ðr À r t ; hÞW 0 ðr À r s ; hÞ r À r s jr À r s j ðm T s m t Þdr : This integration is over a spherical domain (jr À r t j < 2h) around the target particle with radius 2h (the support of the kernel function).Taking into account the symmetric property, the integration result from Eq. ( 36) is a vector parallel to ðr t À r s Þ.
When the two particles are aligned with respect to the magnetic moment, ðm T s m t Þ > 0. For kernel functions, we usually have W ðr; hÞ > 0 and W 0 ðr; hÞ < 0, so the final integration of Eq. ( 36) is a vector pointing in the opposite to the direction of ðr t À r s Þ.
Huang et al. [21] argue that the current loop model violates momentum conservation.Here, we prove that the force of the source particle on a target is aligned with two particle centers, and this would not lead to extra momentum or angular momentum.Moreover, the force between two particles would be attractive if their magnetic moments are aligned in the same direction.On the surface of the fluid domain, the attractive force from the inner part would drag the surface particle back to the fluid domain, thus leading to an inward-pointing magnetic force.As mentioned above, this inward-pointing property enables the integration of the magnetic force into a more advanced SPH framework such as DFSPH and IISPH, which keeps a stable simulation of ferrofluids with large temporal steps.
The kernel function W ðr; hÞ ¼ 0 if jrj > 2h.Now we have two situations to further deal with the integration in Eqs.(30) and (36).When jr s À r t j > 4h, the non-zero part of the kernel function would have no overlap, and in the case jr s À r t j 4h, there are overlaps.
In the first case, jr s À r t j > 4h.M s ¼ m s W ðr À r s ; hÞ ¼ 0 in the spherical domain jr À r t j < 2h (please note, that jr À r s j > 2h).According to Eq. ( 15), r Á H ¼ Àr Á M s ¼ 0. Thus H is harmonic in this domain(i.e., divergence-free and curl-free).Using the mean-value property, we obtain Z F Kelvin dr ¼ m 0 rHðr t À r s ; m s Þm t ; and Thus, In the second case, jr s À r t j 4h, the integral in Eq. ( 29) is first evaluated numerically for a series of relative positions to form look-up curves.Later in the simulation, we use these look-up curves for the fast integral evaluation.The numerical integration is carried out in a coordinate system where we put the source particle at ð0; 0; 0Þ and the target particle at ð0; 0; qhÞ: 0 q 4 is the normalized distance.The rotation matrix is denoted by R, such that m ¼ R m and m ¼ R T m, where the tilde sign "$" indicates that the variable is represented in local coordinates.
The third order tensor Labg below describes the bilinear relationship of the current loop force in Eq. ( 29) with respect to the source and target moments ms and mt .
Labg ðjr t À r s j; hÞ mb s mg t ; a 2 f1; 2; 3g: The parameter a refers to the force component.b and g refer to the source particle and target particle magnetic moment component direction.There are 27 entries in the third order tensor.With numerical integration of Eqs. ( 30) and ( 36), we calculate the forces from the source to the target with unit magnetic moments but with 3 Â 3 ¼ 9 different direction combinations for a series of distances to obtain 27 curves.However, following the work of Huang et al. [21], only 7 of them are non-zero and there are two types of curves for the Kelvin model.In our new force model, the second integration term R F M dr in Eq. ( 36) is added.A necessary condition for this integral to be non-zero is that ð mT s mt Þ should be non-zero, so the terms only take effects when the magnetic moment of two particles are aligned.Thus three combinations L311;322;333 would be modified and lead to three curves, Each curve is fitted by a piece-wise fourth-order polynomial The coefficients of C 1 ; C 2 and C 3 in four intervals, 0 q < 1, 1 q < 2, 2 q < 3, and 3 q < 4 are given in Tables 1,  2, and 3 respectively.Fig. 1 illustrates these different curves.Two curves C 2 and C 3 have negative values when q < 2 resulting in an attractive force.
After calculating the force in local coordinates, we transform it back to global coordinates with f s!t ¼ R fs!t .For the summed magnetic force on a certain particle, we obtain Regarding the force from the external magnetic field, we can easily compute rH ext ðr i Þ and obtain and the total magnetic force given as  Algorithm 2 provides a summary of the magnetic force computation process and our force model fits well into the DFSPH framework in Algorithm 1.

Algorithm 2. Magnetic Force Computation
Input: Positions r Output: Magnetic force f 1: calculate H ext ðr i Þ 2: solve the linear Eq. ( 25) using conjugate gradient 3: update m i based on Eq. ( 24) 4: for i 1 to N do 5: if jr t À r s j 4h then 8: calculate q jr j À r i j=h 9: calculate mj R T m j 10: compute force fj!i using Eq. ( 40) 11: else jr t À r s j > 4h then 14: compute force f j!i using Eq.(39) 15: The force in our new model has an attractive direction which is also illustrated in Fig. 2. The external magnetic field is constantly pointing upwards and used to magnetize the particles in the box.The solution of the magnetic flux density field is shown.The Kelvin model leads to the force pointing outwards, and our current loop model leads to the inward-pointing direction of the force which prevents the levitation artifacts.

RESULTS
In this section, we present different numerical experiments illustrating the effectiveness of our approach.First, a comparison of our current loop model to the Kelvin model is presented.Then the effect of varying surface tension and magnetic field strength is shown.The interaction of the ferrofluid and a cubic magnet is shown in a rotating magnet and in a moving magnet scene.Finally, the direct contact behavior of a rigid magnet and the ferrofluid is presented in a climbing fluid and in a fluid emitter scene.Our improvement for solving the equations describing the magnetic field is shown in the last part.Table 4 provides an overview of the parameters used in our experiments.

Comparing Force Models
As illustrated in the previous section, and in Figs. 1 and 2, the Kelvin force model results in a force pointing outwards with respect to the fluid surface, while the current loop model force points inwards.In the previous work of Huang et al. [21], they use WCSPH, which enables negative pressure.When the particle goes far away from the fluid surface, the negative pressure will drag it back and leads to a stable fluid surface, thus the particle levitation problem is alleviated.However, the negative pressure is known to lead to unnatural particle clustering [26] in density summation based SPH.Adding Kelvin force into DFSPH and IISPH model will lead to noisy particles levitation.Fig. 3 shows the comparison of the dynamic results of the two models at different time steps for DFSPH and IISPH.We can clearly see that, for both DFSPH and IISPH, the Kelvin model shows very noisy particle levitation, while the current loop model manages to maintain a stable fluid surface.
Please note, that in previous work of Huang et al. [21], the time step size is set to be 3 Á 10 À5 s which is greatly   limited by the nature of WCSPH.While in DFSPH and IISPH, the time step size is set to be 5 Á 10 À4 s, which is 16-times larger than the previous one.Hence, our new current loop model enables the fast and stable simulation of ferrofluids.

Varying Surface Tension and Field Strength
We conduct an experiment with varying surface tension and magnetic field strength.While the external magnetic field mainly changes the amplitude of the spikes, the surface tension mainly changes the width or pattern of the fluid's surface.Fig. 4 shows the corresponding numerical result.
Along the horizontal axis, we change the surface tension coefficients.As we can see, the number of spikes gets smaller when we have larger surface tension.Along the vertical axis, we change the magnetic field strength.With larger field strength, we can see an increase of the amplitude.This observation shows that our model can reproduce this property of ferrofluids.

Rotating Magnet
In the previous two experiments, the spikes are generated under a magnetic field which is constant in space and time.The particle sizes are chosen sufficiently small to capture the structure of the characteristic spikes.The susceptibility values is taken from Huang et al. [21].Using the surface tension model [39], we are choose values from 0.25 to 1 for the surface tension coefficients resulting in proper spike shapes.Our SPH framework enables the adaptive time step from 1 Á 10 À4 s to 1 Á 10 À3 s.The runtimes of our numerical experiments are evaluated with respect to a single frame.Please note, that we are exporting 200 frames per second.Here, we present a more complex interaction between the ferrofluid and a rigid cube magnet.When the cube magnet approaches the fluid, it will magnetize the fluid particles.The spikes usually grow in the direction of magnetic field lines.In Fig. 5, the right side has a magnetic field in the vertical direction, such that spikes grow in the vertical direction.While on the left, the magnetic field direction is parallel to the fluid bottom such that the spike will not grow and just cluster together.

Moving Magnet
While the external static magnet can attract the fluid and generate spikes, in this moving magnet experiment, we further show the capability of our framework.In Fig. 6, an external rigid cube magnet is used to generate spikes on the fluid's surface.By slowly moving the external magnet from left to right, the fluid follows the magnet.During the pathway, it interacts with a non-magnetic cylinder.It splits apart when hitting the cylinder and coalesces after leaving it.Our work is able to capture the dynamical process.

Climbing Fluid
In this experiment, we show the direct interaction of the ferrofluid and an external rigid magnet.This hemisphere is magnetized by a cylindrical magnet located below.We first emit fluid over the hemisphere with the magnetic force turned off.After waiting a while with almost resting of the fluid, the magnetic force is turned on.Fig. 7 shows the fluid shape before the magnetic force is turned on and the steady state after the magnetic force turned on.In the dynamic process, we can see that the fluid climbs over the hemisphere and grows spikes upon it.

Fluid Emitter
In this complex scene, a cylinder is carved with helix shape groove.Another cylindrical magnet below is used to magnetize it.Right from the start of the simulation, the ferrofluid is kept been emitted over the helix.The magnetic strengths grows from zero to its maximum at t ¼ 0:5 s.Fig. 8 shows the dynamical result of our method with the rightmost figure showing the state at t ¼ 3:5 s.Compared to other cases above, a smaller time step is used.This is due to the falling down of droplets as their large velocity limits the time step to fulfill the Courant-Friedrichs-Lewy condition.
However, even at the lower boundary, we can still use a time step about three times lager compared to previous    work.When the droplet falls off, it gets attracted by the cylinder and spins along the helix curve down to the bottom.
Our framework is able to capture this complex physical behavior.

Conjugated Gradient
In the work of Huang et al. [21], the magnetic equation is solved by a least squares conjugate gradient solver.However, the linear system of equations for solving the magnetic equation is symmetric.Hence, we can use the standard conjugate gradient solver.Least squares conjugate gradient solves a normal equation, whose condition number is much larger than the original equation, thus it converges slower.In Fig. 9, we present the residual of the two methods.The performance of a standard conjugate gradient is better than we expected.While the least squares conjugate gradient solver could be used in the non-linear magnetization case, standard conjugate gradient turned out to be sufficient in our work when the susceptibility is constant in each scene corresponding to the linear magnetization problem.

CONCLUSION
In our work, we propose a novel magnetic force model.Compared to previous work, our new model avoids the problem of particle levitation, and thus can be nicely incorporated into summation based SPH frameworks, such as IISPH and DFSPH.Compared to WCSPH, our new framework enables a time step about ten times larger.
Although the surface tension model from Akinci et al. [39] enables us to generate proper spike shapes, it does not provide physically accurate surface tension coefficients.To match the exact physical behavior of ferrofluids in a quantitative sense, a more sophisticated surface tension model is required.The Kelvin force model and the current loop force model should theoretically lead to the same surface shape for incompressible ferrofluids, however, their behaviors are slightly different in the particle-based discretized form.The inwardpointing current loop force around the fluid surface tends to smooth the shape of the spikes.
Our framework can also be extended with the fast multipole method for a faster calculation of the magnetic force.However, this is orthogonal to our contrition which is why we did not explore this path in our current work.In our work and Huang et al. [21], SPH is used as the fluid solver.In future work a hybrid Eulerian-Lagrangian fluid solver, such as FLIP, can serve as the backbone of the ferrofluid simulation.The solid-liquid coupling of ferrofluids is also an interesting topic to explore.Fig. 8. From left to right, we show the states of the fluid at a time evolving sequence.Please note, that there is another cylindrical magnet located below which is used to magnetize the helix carved cylinder.Fig. 9.The residual norm between our newly implemented solver compared to the previous one of Huang et al. [21].The blue curve corresponds to the new conjugate gradient solver while the orange curve corresponds to the least squares conjugate gradient solver [21].Here, we are solving for the magnetic field of a sphere with radius of 2 cm.The susceptibility is set to be 500.

Fig. 1 .
Fig. 1.Illustration of the new curves C 1 , C 2 , C 3 in our current loop model, and the old curves C 0 1 , C 0 2 from Huang et al. [21].C M corresponds to the extra term induced by the second integration of Eq. (29).

Fig. 2 .
Fig.2.Sketch of the external magnetic field strength, the magnetic flux density, the corresponding Kelvin force, and the current loop model force on each particle.

Fig. 3 .
Fig. 3.The comparison experiments between the current loop model and Kelvin model are carried out for both DFSPH and IISPH.The external magnetic field strength is set to be a constant value of 1:3 Á 10 5 A/m along the vertical direction and the surface tension coefficient is set to be 0:5.The box size is 0.06 m Â 0.06 m Â 0.015 m.The top two rows refer to the results obtained after t ¼ 0:05 s and the bottom two rows refer to the results obtained after t ¼ 0:10 s.From left to right in each row, the four subplots refer to DFSPH using the current loop model, DFSPH using the Kelvin model, IISPH using the current loop model, and IISPH using the Kelvin model.

Fig. 5 .
Fig.5.Ferrofluids exposed to a rigid cube magnet with a edge length of 4 cm.On the left side, the magnetic field direction is parallel to the fluid bottom, while on the right side the magnetic field direction is vertical pointing up and perpendicular to the fluid bottom.

Fig. 6 .
Fig. 6.The time sequence of this experiment is shown in an N-shape direction.The cube magnet leads to a magnetic field direction along the vertical axis.The edge length of the cube is 4 cm.

Fig. 7 .
Fig. 7. Corresponding rest fluid shapes without (left) and with (right) magnetic force.Please note, that there is another cylindrical magnet located below which is used to magnetize the steel sphere.

TABLE 1 Polynomial
Coefficients of C 1 ðqÞ

TABLE 3 Polynomial
Coefficients of C 3 ðqÞ

TABLE 4 This
Table Provides an Overview of the Experiments Conducted in Our Work