2.1 Set-up and visualization

To define the physical basis for the numerical computations, first we experimentally demonstrate the phenomenon of a non-magnetizable droplet levitation inside a ferrofluid. This approach has actually helped us to first understand about a realistic set of magnetic field boundary conditions that can lead to a stable levitation of the droplet.

Figure 1. (a) The PMMA sheets used to manufacture the Hele-Shaw cell together with a pair of neodymium magnets. (b) Recorded image sequence showing a water droplet levitating inside a ferrofluid sample contained in the Hele-Shaw cell. Only a part of the cell is shown. The droplet radius is ${\approx}2$  mm. The levitation is initiated when the cell is placed on a magnet surface. The gravity acts vertically downward in the figure.

The levitation of the droplet is visualized in a Hele-Shaw cell arrangement. The cell is made up of two closely spaced transparent poly-methyl-methacrylate (PMMA) sheets of size $20~\text{mm}\times 20~\text{mm}$ , shown in figure 1(a). The figure shows the PMMA sheets prior to the manufacture of the cell. The four sides of the cell are closed by inserting a $1.0\pm 0.01~\text{mm}$ thin polymer sheet cuttings in between the PMMA sheets. The final effective volume of the cell is measured to be $17.7~\text{mm}\times 17.7~\text{mm}\times 1.0~\text{mm}$ .

The cell is filled with a sample of kerosene-based ferrofluid (a precise description of the magnetic characterization of the ferrofluid sample, apparatus and image analysis is given later in § 9 while comparing with the simulations). A small droplet of water with volume $\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}_{H}R^{2}$ is placed into the cell before its closure, where $R$ is the droplet radius and $\unicode[STIX]{x1D6FF}_{H}$ is the gap between the walls of the Hele-Shaw cell. The size of the droplet is predicted using image processing after it is dispensed into the cell ( $R\approx 2~\text{mm}$ ). The plane of the cell is kept parallel to the direction of the gravity, and thus in the absence of the field, the droplet remains settled at the bottom.

The levitation against the gravity is now initiated by placing the cell on the surface of a permanent magnet; the magnets used are of neodymium and a pair is shown in figure 1(a) together with the cell. Figure 1 shows the phenomenon recorded during one of the experiments under the above stated conditions. Seven different time instants during the evolution of the droplet are shown in the figure. Initially, the droplet rests on the bottom wall. After the application of the field at the bottom of the cell, the droplet begins to levitate and the shape of the droplet alters – it elongates laterally and attains a concavity at the bottom. Eventually it attains almost an elliptic shape.

2.2 Stable levitation

Under the influence of a single magnet at the base, it is noticed that the location of the levitated droplet is not exactly stable; there is a side-wise movement of the droplet in addition to its vertical rise and the levitation path is not exactly a straight trajectory. In other words, a stable levitation is not achieved using a single magnet at the bottom; after some initial period the droplet path tilts towards the side walls. To achieve a stable levitation and to restrict the side-wise drift of the droplet, an additional pair of neodymium magnets, one at each side wall of the cell, is attached to provide a force on the droplet directing from the side walls towards the centre of the cell. The single magnet at the bottom is now also replaced with a pair so that the north and south ends of the magnets make contact with each other in an alternating fashion and thus remain attached to the cell. That is, if the north pole of a magnet is in contact with the wall, then its adjacent magnets will have their south poles in contact with the respective walls – an arrangement very similar to the Halbach array of magnets (Halbach Reference Halbach1985). One pair of magnets at the top wall closes the loop. The presence of this complimentary pair at the top wall is expected to reduce the levitation height of the droplet, but its absence gives rise to a field configuration which affects the stability of the levitation negatively and in an interesting way (here briefly, the absence of this pair of magnets at the upper wall gives rise to more than two possible equilibrium locations for the droplet; we elaborate on this in § 8). In figure 1(a), one such pair of magnets is shown. The complete schematic for the magnet arrangement is shown in figure 3, and is further explained below in § 3.

Figure 2. The recording of one demonstrative experiment with a weakly magnetizable ferrofluid sample. The bottom pair of magnets remains fixed at the base and the cell is placed onto the magnets, initiating the levitation of the non-magnetizable droplet.

One of the experimental demonstrations conducted with this arrangement of magnets is shown in figure 2. Eight different time instants are shown in the figure. Initially the water droplet remains in the vicinity of the bottom wall. Under this arrangement of magnets, the droplet eventually finds a stable equilibrium position depending upon the balance between the gravitational, buoyancy and the magnetic forces. Certainly, other magnet arrangements for a stable levitation of the water droplet may also be possible; the current alternating arrangement has proven specifically helpful in incorporating the magnetic field boundary conditions conveniently in terms of sine functions multiplied by the maximum strength of the magnet (§ 3.2). This Halbach array arrangement is also practically helpful in keeping the magnet edges in contact with each other without any external forcing.

2.3 Physical explanation

A possible physical explanation of the phenomena can be given after Rosensweig (Reference Rosensweig1985). First, let us assume that there is no droplet. When the cell is brought into contact with the magnet, and then held stationary, the ferrofluid is attracted towards the magnet. For an initial period of time, some transients appear inside the ferrofluid due to this attractive force. However, as the cell is then kept stationary and there is not a continuous supply of work, the stationary magnetic force cannot make the ferrofluid move continuously due to the thermodynamic constraints. The stresses inside the ferrofluid reorient themselves to counter this perpetual motion, much like the pressure redistribution in the gravitational field. The net response from the ferrofluid in this static condition is the redistribution of the mechanical pressure to balance the magnetic and gravitational forces; the pressure being higher in higher magnetic field regions. Now, if a non-magnetizable object is placed inside the ferrofluid, it experiences the developed pressure gradient and starts moving away from the magnet. In the next section we describe the mathematical basis of our computations.

3.1 Dimensional form of the governing equations

The motion of the non-magnetic droplet inside the ferrofluid medium, and the continuous local flow and the magnetic field, are described by a coupling between the Navier–Stokes equations and Maxwell’s equations of electromagnetism. If the fluids are considered electrically non-conducting, and assuming that the relaxation time of the magnetic nano-particles in the ferrofluid is much smaller than the relevant hydrodynamic time scales, then Maxwell’s equations reduce to the magnetostatic form,

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0,\quad \unicode[STIX]{x1D735}\times \boldsymbol{H}=0,\\ \boldsymbol{B}=\unicode[STIX]{x1D707}_{o}(\boldsymbol{M}+\boldsymbol{H}),\end{array}\right\}\quad \text{in}~\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\boldsymbol{B}$ , $\boldsymbol{H}$ and $\boldsymbol{M}$ are the magnetic flux density, the magnetic field and the magnetization respectively while $\unicode[STIX]{x1D707}_{o}$ is the free space magnetic permeability. The irrotationality of $\boldsymbol{H}$ permits the relation

(3.2) $$\begin{eqnarray}\boldsymbol{H}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is a scalar magnetic potential. Also $\boldsymbol{M}$ is constitutively related to $\boldsymbol{H}$ through

(3.3) $$\begin{eqnarray}\boldsymbol{M}=\unicode[STIX]{x1D712}\boldsymbol{H},\end{eqnarray}$$

where $\unicode[STIX]{x1D712}$ is the magnetic susceptibility. For ferrofluids in general, the magnetization is itself governed by a differential equation involving the magnetization relaxation time, however, here we assume that the relaxation time is small and the magnetization relaxes in an infinitesimally small time, the so-called quasi-equilibrium ferrohydrodynamic hypothesis. A single scalar equation can be obtained using (3.2) and (3.3) in (3.1), expressed as

(3.4) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}_{o}(1+\unicode[STIX]{x1D712})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0,\quad \text{in}~\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where the quantity $\unicode[STIX]{x1D707}_{o}(1+\unicode[STIX]{x1D712})$ is equal to the magnetic permeability $\unicode[STIX]{x1D707}$ .

The ferrofluid susceptibility varies with the magnetic field strength and other thermodynamic variables. In this study we consider it either a constant or a function of the local magnetic field. The simplifying assumption of constant permeability of the ferrofluid is used in the first two sets of simulations to obtain the basic characteristics of the droplet shape and levitation height, while being computationally efficient. The results from the constant permeability assumption have also served as a reference for further refinement of our simulations using a more realistic field-dependent permeability model, when the experimentally observed shape is compared with the simulations. In the latter case, the magnetization relation (3.3) can be written in its nonlinear form, i.e. $\boldsymbol{M}=\unicode[STIX]{x1D712}(H)\boldsymbol{H}$ . Then the equation for the magnetic potential is different for $\unicode[STIX]{x1D6FA}_{f}$ and $\unicode[STIX]{x1D6FA}_{d}$ and is expressed as (using Langevin’s nonlinear magnetization equation)

(3.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}_{o}\left(1+\frac{M_{s}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\left[\coth \unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|-\frac{1}{\unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\right]\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0,\quad \text{in}~\unicode[STIX]{x1D6FA}_{f}~\text{if}~\unicode[STIX]{x1D712}_{f}=\unicode[STIX]{x1D712}_{f}(H),\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}_{o}(1+\unicode[STIX]{x1D712}_{f})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0,\quad \text{in}~\unicode[STIX]{x1D6FA}_{f}~\text{if}~\unicode[STIX]{x1D712}_{f}~\text{is constant},\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}_{o}(1+\unicode[STIX]{x1D712}_{d})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0,\quad \text{in}~\unicode[STIX]{x1D6FA}_{d},\end{array}\right\}\end{eqnarray}$$

where $M_{s}$ is the saturation magnetization of the ferrofluid and $\unicode[STIX]{x1D6FE}=3\unicode[STIX]{x1D712}_{o}/M_{s}$ . It should be noted that this specific nonlinear ferrofluid magnetization model assumes that the ferrofluid exhibits nearly a paramagnetic behaviour and also does not reflect the hysteresis of magnetization. It further assumes a nearly mono-disperse size distribution of the magnetic nano-particles and negligible dipole–dipole interactions. For non-magnetic liquids, the magnetic susceptibility is usually negligible, or in other words, the permeability can be considered equal to the free space permeability  $\unicode[STIX]{x1D707}_{o}$ . Equation (3.5) serves as the key to obtaining $\boldsymbol{H}$ over the whole domain  $\unicode[STIX]{x1D6FA}$ .

The equations for the isothermal and incompressible flow field are

(3.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0,\\ \displaystyle \unicode[STIX]{x1D70C}\frac{\boldsymbol{D}\boldsymbol{v}}{\boldsymbol{D}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D64E}+\unicode[STIX]{x1D64E}_{m})+\unicode[STIX]{x1D70C}\boldsymbol{g}+\boldsymbol{f}_{s},\end{array}\right\}\quad \text{in}~\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $p,\unicode[STIX]{x1D64E},\unicode[STIX]{x1D64E}_{m},\boldsymbol{g}$ and $\boldsymbol{f}_{s}$ are the mechanical pressure, viscous stress tensor, magnetic stress tensor, gravitational acceleration and the interfacial force respectively. The Newtonian viscous stress tensor is $\unicode[STIX]{x1D702}(\unicode[STIX]{x1D735}\boldsymbol{v}+(\unicode[STIX]{x1D735}\boldsymbol{v})^{\text{T}})$ . The force $\boldsymbol{f}_{s}$ is expressed as $\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}\unicode[STIX]{x1D6FF}_{s}$ and acts singularly at the interface. Here $\unicode[STIX]{x1D705}$ is the local curvature of the interface, $\boldsymbol{n}$ is the unit outward normal at the interface and $\unicode[STIX]{x1D6FF}_{s}$ is the delta function at the interface.

The driving force in the present study is due to the magnetic stresses. There are number of different expressions for $\unicode[STIX]{x1D64E}_{m}$ that exist in the ferrohydrodynamic literature; interestingly, under incompressible and isothermal conditions and for isotropic permeability $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}(H)$ , the existing expressions for the magnetic stress tensor reduce to the form

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{m}=-a\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}\boldsymbol{H}\boldsymbol{H}^{\text{T}},\end{eqnarray}$$

where the first part $-a\unicode[STIX]{x1D644}$ , which contains the isotropic magnetic pressure, can safely be lumped with $-p\unicode[STIX]{x1D644}$ (Rosensweig Reference Rosensweig1985; Afkhami et al. Reference Afkhami, Renardy, Renardy, Riffle and St Pierre2008). This form of the magnetic stress tensor is particularly well adaptable to a conservative finite-volume formulation in comparison to the magnetic body force density expressions such as the Kelvin force density or the Korteweg–Helmholtz force density which are derivable by taking the divergence of  $\unicode[STIX]{x1D64E}_{m}$ . We directly discretize the divergence of (3.7) on finite-volume cells utilizing $\int _{V}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{m}\,\text{d}V=\int _{S}\boldsymbol{n}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{m}\,\text{d}S$ (where $V,S$ and $\boldsymbol{n}_{s}$ denote the computation cell volume, cell surface and the outward normal at the cell surface respectively), which conserves the magnetic force fluxes for each individual computational cell.

3.2 Initial, boundary and interfacial conditions

The no-slip and the no-penetration flow boundary conditions are considered at all four walls while sinusoidal boundary conditions are considered for the gradient of magnetic potential to replicate the Halbach array of magnets, described as

(3.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{v}=\mathbf{0},\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA},\\ \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\boldsymbol{n}_{b}=\left\{\begin{array}{@{}l@{}}H_{o}\sin (2\unicode[STIX]{x03C0}x/L),\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{t},\\ H_{o}\sin (2\unicode[STIX]{x03C0}y/L),\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{l},\\ -H_{o}\sin (2\unicode[STIX]{x03C0}y/L),\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{r},\\ -H_{o}\sin (2\unicode[STIX]{x03C0}x/L),\quad \text{on}~\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{b},\end{array}\right.\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{n}_{b}$ is the unit normal at the walls pointing into $\unicode[STIX]{x1D6FA}$ and subscripts $t,l,r$ and $b$ represent top, left, right and bottom wall respectively.

Although the conditions at the interface are not explicitly needed in the numerical treatment of the two-phase, one-fluid formulation (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011), we state them below for the sake of completeness. At the interface, the normal component of the magnetic flux density, the tangential component of the magnetic field, the tangential component of the total stress and the velocity are continuous, while there is discontinuity in the normal component of the total stress,

(3.9) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{n}\boldsymbol{\cdot }[\boldsymbol{B}]=0,\quad \boldsymbol{n}\times [\boldsymbol{H}]=0,\quad [\boldsymbol{t}^{\text{T}}\unicode[STIX]{x1D64F}\boldsymbol{n}]=0,\quad [\boldsymbol{v}]=0,\\ \text{}[\boldsymbol{n}^{\text{T}}\unicode[STIX]{x1D64F}\boldsymbol{n}]=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D64E}+\unicode[STIX]{x1D64E}_{m}$ is the total stress tensor and $[x]$ denotes the difference of a quantity, $x$ , right across the interface.

3.3 Non-dimensionalization

From the governing model, the flow solution is dependent on the following dimensional parameters

(3.10) $$\begin{eqnarray}\boldsymbol{v}=\boldsymbol{v}(\boldsymbol{x},t;\unicode[STIX]{x1D70C},\unicode[STIX]{x1D702},R,g,\unicode[STIX]{x1D70E},\unicode[STIX]{x1D707},\unicode[STIX]{x1D6FE},H_{o},M_{s}).\end{eqnarray}$$

The number of independent variables is reduced by properly identifying the non-dimensional groups which influence the flow solution. Taking $\unicode[STIX]{x1D6FA}_{d}$ as the reference for all of the properties except the magnetic permeability, for which ferrofluid has a permeability higher than that of the droplet medium, and considering the following reference scales,

(3.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\sim R^{-1},\quad \unicode[STIX]{x1D70C}\sim \unicode[STIX]{x1D70C}_{d},\quad \unicode[STIX]{x1D702}\sim \unicode[STIX]{x1D702}_{d},\quad \unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{f},\\ \boldsymbol{v}\sim \unicode[STIX]{x1D702}_{d}/\unicode[STIX]{x1D70C}_{d}R,\quad t\sim \unicode[STIX]{x1D70C}_{d}R^{2}/\unicode[STIX]{x1D702}_{d},\quad p\sim \unicode[STIX]{x1D702}_{d}^{2}/\unicode[STIX]{x1D70C}_{d}R^{2},\\ \boldsymbol{H}\sim H_{o},\quad \boldsymbol{B}\sim \unicode[STIX]{x1D707}_{o}H_{o},\quad \boldsymbol{M}\sim M_{s},\quad \boldsymbol{f}_{s}\sim \unicode[STIX]{x1D70E}/R^{2},\end{array}\right\}\end{eqnarray}$$

the equations of fluid motion are normalized to the following form

(3.12) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\boldsymbol{v}^{\ast }=0,\\ \displaystyle \unicode[STIX]{x1D70C}^{\ast }\frac{\boldsymbol{D}\boldsymbol{v}^{\ast }}{\boldsymbol{D}t^{\ast }}=\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }(-p^{\ast }\unicode[STIX]{x1D644}+\unicode[STIX]{x1D64E}^{\ast })+La_{m}\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{m}^{\ast }+Ga\,\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{g}^{\ast }+La\,\boldsymbol{f}_{s}^{\ast },\end{array}\right\}\quad \text{in}~\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where the non-dimensional group

(3.13) $$\begin{eqnarray}La=\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D70C}_{d}R}{\unicode[STIX]{x1D702}_{d}^{2}}\end{eqnarray}$$

is the Laplace number signifying the ratio of the interfacial force to the viscous force,

(3.14) $$\begin{eqnarray}La_{m}=\frac{\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D70C}_{d}H_{o}^{2}R^{2}}{2\unicode[STIX]{x1D702}_{d}^{2}}\end{eqnarray}$$

is the magnetic Laplace number signifying the ratio of the magnetic force to the viscous force and

(3.15) $$\begin{eqnarray}Ga=\frac{g\unicode[STIX]{x1D70C}_{d}^{2}R^{3}}{\unicode[STIX]{x1D702}_{d}^{2}}\end{eqnarray}$$

is the Galilei number signifying the ratio of the gravitational force to the viscous force.

Maxwell’s equations transform to the following non-dimensional form,

(3.16a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\boldsymbol{B}^{\ast }=0,\quad \unicode[STIX]{x1D735}^{\ast }\times \boldsymbol{H}^{\ast }=0,\quad \boldsymbol{B}^{\ast }=\boldsymbol{H}^{\ast }+\left[\frac{1}{\unicode[STIX]{x1D709}_{o}}\right]\boldsymbol{M}^{\ast },\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{o}=H_{o}/M_{s}$ . The magnetization relation takes the following normalized form

(3.17) $$\begin{eqnarray}\boldsymbol{M}^{\ast }(\boldsymbol{H}^{\ast })=\displaystyle \left[\coth \unicode[STIX]{x1D6FE}_{o}H^{\ast }-\frac{1}{\unicode[STIX]{x1D6FE}_{o}H^{\ast }}\right]\frac{\boldsymbol{H}^{\ast }}{H^{\ast }}=\unicode[STIX]{x1D712}^{\ast }(\unicode[STIX]{x1D6FE}_{o},H^{\ast })\boldsymbol{H}^{\ast },\end{eqnarray}$$

where

(3.18) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{o}=\frac{3\unicode[STIX]{x1D712}_{o}H_{o}}{Ms}=3\unicode[STIX]{x1D712}_{o}\unicode[STIX]{x1D709}_{o}.\end{eqnarray}$$

Using this relation for magnetization, the equation for magnetic potential takes the following non-dimensional form

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|}\left[\coth \unicode[STIX]{x1D6FE}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|-\frac{1}{\unicode[STIX]{x1D6FE}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|}\right]\right)\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }=0,\end{eqnarray}$$

or

(3.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}}\unicode[STIX]{x1D712}^{\ast }(\unicode[STIX]{x1D6FE}_{o},H^{\ast })\right)\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }=0.\end{eqnarray}$$

Therefore, four non-dimensional parameters, $La,La_{m},Ga$ and $\unicode[STIX]{x1D6FE}_{o}$ are significant and thus the flow solution in the case of field-dependent ferrofluid susceptibility can be expressed by the functional form

(3.21) $$\begin{eqnarray}\boldsymbol{v}^{\ast }=\boldsymbol{v}^{\ast }(\boldsymbol{x}^{\ast },t^{\ast };La,La_{m},Ga,\unicode[STIX]{x1D6FE}_{o}).\end{eqnarray}$$

In the case of constant ferrofluid susceptibility, where there is no bound on  $\boldsymbol{M}$ , the reference scale for the magnetization is chosen equal to  $H_{o}$ . In such a case the functional form of the flow solution reads

(3.22) $$\begin{eqnarray}\boldsymbol{v}^{\ast }=\boldsymbol{v}^{\ast }(\boldsymbol{x}^{\ast },t^{\ast };La,La_{m},Ga).\end{eqnarray}$$

The parameter $\unicode[STIX]{x1D6FE}_{o}$ depicts the effect of nonlinearity of the magnetization on the flow solution. In the following sections we drop the star symbols on the non-dimensional variables for convenience.

4.1 Description of the algorithm

Initially at $t=0$ , both the fluids are considered quiescent $(\boldsymbol{v}(\boldsymbol{x},0)=\mathbf{0})$ . The initial location of the interface between the two phases is considered to be known, or in other words, the initial discrete distribution of the fluid property fields – $\unicode[STIX]{x1D70C}_{i,j},\unicode[STIX]{x1D702}_{i,j},\unicode[STIX]{x1D707}_{i,j}$ – on the grid is considered to be known; $i,j$ being the indices associated with the concerned grid point. Knowing $\unicode[STIX]{x1D707}_{i,j}$ at $t=0$ , the equation for the magnetic potential $\unicode[STIX]{x1D719}$ (3.4) is first solved. Obtaining $\unicode[STIX]{x1D719}$ , $\boldsymbol{H}$ is obtained from $\boldsymbol{H}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . Once $\boldsymbol{H}$ is known, the magnetic body force is computed. The equation of motion (3.6) is then solved for the velocity field $\boldsymbol{v}$ using the pressure projection algorithm (Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011). Using $\boldsymbol{v}$ , the location of the interface is then advected, and from this updated location of the interface, the discrete marker function $C_{i,j}$ is reconstructed using the front-tracking algorithm (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992; Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011). The discrete density and viscosity fields are then interpolated from $C_{i,j}$ as

(4.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{i,j}=\unicode[STIX]{x1D70C}_{f}C_{i,j}+\unicode[STIX]{x1D70C}_{d}(1-C_{i,j}),\quad \unicode[STIX]{x1D702}_{i,j}=\unicode[STIX]{x1D702}_{f}C_{i,j}+\unicode[STIX]{x1D702}_{d}(1-C_{i,j}).\end{eqnarray}$$

The permeability $\unicode[STIX]{x1D707}_{i,j}$ is obtained in a similar fashion if it is assumed constant in both the phases. However, when the nonlinear magnetization model is considered, the permeability is considered constant only in the non-magnetizable phase. In the magnetizable phase, it is predicted using Langevin’s relation for the ferrofluid magnetization. Both of the above cases are summarized as

(4.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D707}_{i,j}=\unicode[STIX]{x1D707}_{f}C_{i,j}+\unicode[STIX]{x1D707}_{d}(1-C_{i,j}),\quad \text{if}~\unicode[STIX]{x1D707}_{f}~\text{is assumed constant},\\ \displaystyle \unicode[STIX]{x1D707}_{i,j}=\unicode[STIX]{x1D707}_{o}\left(1+\frac{M_{s}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\left[\coth \unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|-\frac{1}{\unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\right]\right)C_{i,j}+\unicode[STIX]{x1D707}_{d}(1-C_{i,j}),\quad \text{if}~\unicode[STIX]{x1D707}_{f}=\unicode[STIX]{x1D707}_{f}(H).\end{array}\right\}\end{eqnarray}$$

Once the property fields are advected, the algorithm is advanced to the next time step and the numerical cycle is repeated until the desired time.

The use of constant magnetic permeability is a greater computational simplification and results in lower computational times. Therefore it is useful for obtaining very basic features of the flow, and also when parameters other than magnetic permeability (or susceptibility) are varied. An example is the variation of the viscosity ratio, which is not related to the permeability or magnetization of the ferrofluid in our formulation. However, instead of varying the permeability ratio, and for comparison with experiments, we adopt the nonlinear magnetization model.

The spatial as well as the temporal discretizations in our scheme are second-order accurate. The advection term in the momentum equation is handled using a second-order essentially non-oscillatory (ENO) scheme while the standard second-order centred in space discretization is applied to the diffusive viscous terms. The interpolations of the viscosity at the computational cell faces are performed harmonically while the density and the permeability are interpolated arithmetically. A highly optimized V-cycle MULTIGRID method is implemented for the solution of the pressure Poisson equation and the equation for the magnetic potential.

4.2 Set-up of simulations

The properties of the two fluids are presented in table 1, where $\unicode[STIX]{x1D70C}_{d}\sim 10^{3}~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D702}_{f}\sim 10^{-2}~\text{Pa}~\text{s}$ and $\unicode[STIX]{x1D70E}\sim 10^{-3}~\text{N}~\text{m}^{-1}$ . We consider a water droplet of radius of ${\sim}1~\text{mm}$ . The magnetic fields used in our experiments are of $H_{o}\sim 10^{2}~\text{kA}~\text{m}^{-1}$ while the order of permeability of the ferrofluid is considered a multiple of  $\unicode[STIX]{x1D707}_{o}$ . For these orders of magnitude of the properties, the $La\sim 10$ and $Ga$ can be as high as ${\sim}10^{2}$ , while a realistic $La_{m}$ turns out to be ${\sim}10^{4}$ . Considering these, we simulate for $La$ and $Ga$ in the range 0.1–100 while for $La_{m}$ in the range 10–1000. Being conservative, the $La_{m}$ is not increased further to avoid any numerical instability arising due to the higher order of the magnetic source term in the momentum equation.

Table 1. Physical properties of both the phases. The saturation magnetization of the ferrofluid sample is measured by an EverCool SQUID VSM DC magnetometer and the interfacial tension is determined by a technique used in Zhu et al. (Reference Zhu, Nguyen, Ramanujan and Huang2011a ).

a Computed using Langevin’s function for ferrofluid magnetization.

The value of the parameter $\unicode[STIX]{x1D6FE}_{o}$ changes due to changes in  $H_{o}$ , as $\unicode[STIX]{x1D712}_{o}$ and $M_{s}$ do not vary for a given ferrofluid sample. Here, the initial susceptibility of ferrofluids is measured to be ${\approx}0.0819$ , while $M_{s}\sim 10^{3}~\text{A}~\text{m}^{-1}(57.7~\text{G})$ . Thus for an applied field in the range $10^{3}-10^{4}~\text{A}~\text{m}^{-1}$ , the order of $\unicode[STIX]{x1D6FE}_{o}=3\unicode[STIX]{x1D712}_{o}H_{o}/M_{s}$ is expected to vary in the range 0.2–2.0. We simulate for $\unicode[STIX]{x1D6FE}_{o}$ in the range 0.1–2.

The breakup of the simulations and all the corresponding simulation parameters are summarized in table 2. The simulations are divided into six different sets with different foci. In the first set, we simulate for constant permeability in both the phases. Next, the effect of viscosity ratio is investigated. In the third set of simulations, the effects of field-dependent ferrofluid permeability, or equivalently the effects of nonlinear magnetization, are simulated. The stability of levitation is the focus of the subsequent set, while the appearance of interfacial singularity and comparison with experiment are addressed through the last set of simulations. The same hierarchy, as listed in table 2, is utilized for sectioning the results.

Table 2. Breakup of simulation results and corresponding values of input parameters.

a $\unicode[STIX]{x1D6FE}_{o}=3\unicode[STIX]{x1D712}_{o}H_{o}/M_{s}=3\unicode[STIX]{x1D712}_{o}\unicode[STIX]{x1D709}_{o}$ .

b The vacuum permeability is used to compute characteristic $La_{m}$ , as the permeability of the ferrofluid in this case is a variable determined by Langevin’s function.

c The quantity is computed using Langevin’s function for the field-dependent ferrofluid magnetization.

d The quantity is not uniquely defined, or does not appear in the formulation.

4.3 Output fields and variables

Besides the magnetic and velocity field solutions, we study the shapes of the levitating droplet, the final levitation height, the time-dependent displacement and velocity of the droplet tip/nose, the droplet deformation (defined below) and the time-averaged change in the droplet deformation. In certain cases, we also utilize the absolute magnetic field contours and the vorticity contours.

The discrete information about the interface is described by two vectors – $\{x_{k}\}=(x_{1},x_{2},\ldots ,x_{N})^{\text{T}}$ and $\{y_{k}\}=(y_{1},y_{2},\ldots ,y_{N})^{\text{T}}$ – storing $x$ and $y$ coordinates of $N$ discrete interface points. This structure is compatible with the implementation of the front-tracking algorithm. The total number of interface points  $N$ , however, can be different for each time step due to their dynamic addition or deletion, a part of the front-tracking algorithm. From this information about the interface, the vertical distance of the tip of the droplet from the bottom wall is extracted using

(4.4) $$\begin{eqnarray}h_{tip}(t)=\max (\{y_{k}\}),\end{eqnarray}$$

and the speed of the tip is computed using

(4.5) $$\begin{eqnarray}v_{tip}(t)=[h_{tip}(t+\unicode[STIX]{x0394}t)-h_{tip}(t)]/\unicode[STIX]{x0394}t,\end{eqnarray}$$

where $\unicode[STIX]{x0394}t$ is the time step.

To study the time-dependent deformation of the levitating droplet, a global droplet deformation parameter or shape factor, ${\mathcal{D}}$ , is computed from the interface information using

(4.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{D}}(t)=\frac{A_{x}(t)-A_{y}(t)}{A_{x}(t)+A_{y}(t)},\\ \displaystyle A_{y}(t)=\max (\{y_{k}\})-\min (\{y_{k}\}),\\ \displaystyle A_{x}(t)=\max (\{x_{k}\})-\min (\{x_{k}\}).\end{array}\right\}\end{eqnarray}$$

Notice that ${\mathcal{D}}(t)$ can be negative if $A_{x}(t)<A_{y}(t)$ , which implies that the vertical span of the droplet is greater than its horizontal span.

To show the complete ${\mathcal{D}}(t)$ curve for all the simulations is a cumbersome task. To mark the average extent of the droplet deformation over all times, a time-averaged shape factor, $\langle {\mathcal{D}}\rangle$ , is computed. It reduces the results for the shape factor to one value for one simulation. The expression for $\langle {\mathcal{D}}\rangle$ is

(4.7) $$\begin{eqnarray}\displaystyle \langle {\mathcal{D}}\rangle =\frac{1}{t_{end}}\mathop{\sum }_{t=0}^{t_{end}}({\mathcal{D}}(t+\unicode[STIX]{x0394}t)-{\mathcal{D}}(t)),\end{eqnarray}$$

where $t_{end}$ is the time up to which the flow is simulated.

4.4 Code validation and grid and time-step independence

The computational code FERRO is developed and tested for single-phase as well as two-phase viscous incompressible flow, with and without the magnetic effects, and has been successfully applied to problems by the authors to predict the interfacial (Singh, Das & Das Reference Singh, Das and Das2016b ) and relaxation mechanisms (Singh, Das & Das Reference Singh, Das and Das2016a ) in ferrofluids. We incorporate the nonlinear magnetization/field-dependent permeability model and the conservative discretization of the full magnetic body stress tensor, and a highly optimized multigrid technique for the solution of the pressure Poisson and the magnetic potential equations. For the problem at hand, the grid as well as the time step independence of the code are carefully studied. The details of these consistency checks have been presented in the Appendix. Besides these quantitative consistency checks, the physical correctness is also confirmed in our study through its ability to predict the experiments, especially the fine features at the interface (§ 9).

5.1 The effect of change in $La$

The non-dimensional number $La$ is varied for fixed $Ga$ and $La_{m}$ . The interface of the levitating droplet for $Ga=1.0$ and $La_{m}=360$ is shown in figure 5. The flow is simulated for three distinct $La$ – 0.1, 1.0 and 10.0. Note that only a part of the domain ( $0.325\leqslant x\leqslant 0.675$ , $0.0\leqslant y\leqslant 1.0$ ) is shown. The absolute field strength contours are also presented. Whereas the steady state levitation height of the droplet seems nearly independent of $La$ , the difference lies in the extent of the deformation of the interface. For $La=10.0$ , it is only during the very start of the phenomenon (instant $t=0.5$ and 1.0) when the droplet shape deviates from its initial round configuration. The droplet regains a nearly round shape at around $t=1.5$ . As $La$ is decreased, the deformation of the interface becomes more prominent. The tail of the levitating droplet deforms to have a concavity. As time progresses, this feature at the bottom of the droplet disappears in the case of $La=1.0$ . However, this feature at the tail of the droplet stays in the case of $La=0.1$ . The shape in the case of $La=0.1$ seems to resemble a crescent, except for the fact that round projections appear instead of sharp cusps at the tail. The increase in the droplet deformation due to a decrease in $La$ is expected, as a smaller $La$ signifies a lower interfacial energy. However, the transition specifically towards the crescent-like shapes is worth noting. The droplet retains its symmetry about the vertical axis and is asymmetric about the horizontal axis. Its shape is nearly an oval in the case of $La=10.0$ and at low time $t=0.5$ . Lowering the interfacial energy further by reducing $La=1.0$ results in a transition from oval towards crescent. At longer time, the droplet regains its round shape at high $La$ while the deformed shapes persist for low $La$ . The change in the final average levitation height of the droplet due to the change in $La$ is meagre.

Figure 6. The interface of the levitating droplet for different $La_{m}$ . $Ga=1.0$ , $La=0.1$ .

5.2 The effect of change in $La_{m}$

The deformation of the droplet is not solely sensitive to $La$ . In this subset of simulations, the $La$ and $Ga$ numbers are fixed while the deformation is studied with changing  $La_{m}$ . These results are depicted in figure 6. The droplet rises to a relatively greater height and its deformation is also enhanced at increased  $La_{m}$ . A further difference is noted for $La_{m}=640$ and 1000. In these two cases, the shape of the rising droplet is more skirted for $t\leqslant 2.0$ . For $La_{m}=1000$ , the interfacial deformation is more pronounced and the droplet eventually gains a segmented-ring-like shape for $t\geqslant 1.0$ . The increase of $La_{m}$ from 40 to 1000 leads to a further shape transition at late times – first from oval to nearly crescent, and then to a nearly segmented ring.

The levitation height increases with $La_{m}$ due to the stronger field gradients caused by the increase in  $La_{m}$ . The reasons behind the transitions in the shape of the droplet, however, seem more involved. For $La_{m}=1000$ , the initial transition from the round towards a skirted configuration, and eventually to the segmented-ring shape, is particularly intriguing. A possible physical explanation is as follows. In figure 6, the droplet at $t=5.0$ , for all  $La_{m}$ , is close to its steady state. It is clear that as $La_{m}$ increases, the droplet attains its equilibrium increasingly close to the centre of the domain ( $y=0.5$ ), where the magnetic field lines (not absolute field contours) are more distorted. In the case of lower interfacial energy for lower  $La$ , the field patterns near the centre of the domain tend to deform the droplet into a segmented-ring-like shape. On the other hand, lower $La_{m}$ magnitudes fail to raise the droplet to such field line regions, and this results in less alteration in the interface configuration. The magnetic field is primarily responsible for the deformation of the droplet at steady state when the flow velocity vanishes. The hydrodynamic flow field is thus at play primarily during the rise of the droplet before the steady state. As the maximum extent of the concavity at the tail of the droplet is observed during the rise of the droplet, the hydrodynamic flow helps in enhancing the deformation caused by the magnetic interfacial force. Once the steady state is reached, the flow field dies out and only the balance between magnetic, interfacial tension and gravity forces governs the deformed shape. This insight is also quantitatively supported by the time curves of the droplet deformation $\langle {\mathcal{D}}\rangle$ (see the Appendix), where a peak in $\langle {\mathcal{D}}\rangle$ is observed before the steady state. The deformation relatively decreases after this hydrodynamically dominated regime.

The extent of deformation of the droplet with respect to $La$ and $La_{m}$ is quantified in figure 7 with the help of the time-averaged droplet deformation parameter. This supports the fact that the deformation significantly increases with increasing  $La_{m}$ , provided $La$ is low enough. For example at $La=0.1$ , the deformation parameter increases from 0.1 to 0.33 as $La_{m}$ is increased from 40 to 1000. On the other hand, at three orders of magnitude higher $La=10$ , the value of the deformation parameter is negligible, and stays close to zero even for high  $La_{m}$ . This clearly suggest that the deformation of the droplet is maximized at low $La$ and higher  $La_{m}$ . In addition, the time-averaged deformation increases nearly linearly with  $La_{m}$ . The slope of the linear dependence is dependent on  $La$ . At high  $La$ , the slope is close to 0. A decrease in $La$ amplifies the effect of $La_{m}$ on the global deformation of the droplet.

Figure 7. The change in the global time-averaged deformation parameter $\langle {\mathcal{D}}\rangle$ (equation (4.7)) with respect to $La_{m}$ and  $La$ , quantitatively signifying that the deformation of the droplet increases with increasing $La_{m}$ but decreasing  $La$ . $Ga=1.0$ .

5.3 The effect of change in $Ga$

The $Ga$ number signifies the extent of the gravitational force on the droplet relative to the viscous force. We vary $Ga$ for fixed $La$ and $La_{m}$ . Figure 8 depicts the results when $Ga$ is varied from 0.1 to 10.0 at intermediate values of $La=1.0$ and $La_{m}=360$ . As expected, increased $Ga$ decreases the levitation height of the droplet. In addition, the droplet is more suppressed in the case of $Ga=10.0$ as compared to $Ga=0.1$ and 1.0. Observing this, it is certainly possible that if $Ga$ is increased further, the droplet may settle down instead of levitating. Although there is not a considerable deviation of the droplet shape for $Ga=0.1$ and 1.0 from the previously discussed simulations, a new transition to a tooth-like shape occurs near the steady state in the case of higher $Ga$ ( $=10$ ). In this case, if the top portion of the droplet interface becomes concave upward, the shape will nearly be a tooth surface. We show in § 5.5 that for certain combinations of $La$ , $La_{m}$ and $Ga$ , the droplet closely resembles such a shape.

Figure 8. The interface of the levitating droplet for different $Ga$ . $La=1.0$ , $La_{m}=360$ .

5.4 Levitation height and speed characteristics

In this subsection, we study the levitation height and the speed of the droplet as a function of time and look for different temporal modes to reach the equilibrium shape. For the time-dependent levitation height, the vertical location of the tip/nose of the droplet, denoted as $h_{tip}(t)$ (§ 4.3), is tracked with time. The results are depicted in figure 9. The $La$ number increases from left to right while the $Ga$ number increases from top to bottom. In each of the nine plots, five different $La_{m}$ cases are considered.

Figure 9. The time-dependent vertical location of the tip of the droplet for different $La,Ga$ and $La_{m}$ .

It is observed that the $La$ number has the least effect on the $h_{tip}(t)$ curves. Although the droplet shape has turned out to be quite sensitive to $La$ (as discussed previously), its evolution in time shows an opposite nature.

The effect of $La_{m}$ on the $h_{tip}(t)$ curves is quite intriguing. First, we take the case when both $La$ and $Ga$ are 0.01 (top-left plot in figure 9). For $La_{m}$ values of 40 and 160, the droplet approaches an equilibrium location, without crossing the steady state levitation height. However, with further increase of $La_{m}$ to 360 and 640, there start to appear a cross-over/overshoot in the curve. At $La_{m}=1000$ , the overshoot is quite apparent; the value of $h_{tip}$ starts increasing, reaches a maximum and then begin to decrease while approaching the stationary state. Physically this implies that if the order of magnitude of the magnetic levitation force on the non-magnetic droplet, relative to the viscous and gravitational forces on it, is considerably higher, then the droplet can cross the equilibrium location before finally attaining it. This also suggest that the effect of the hydrodynamic flow during the rise of the droplet is not always the same; the viscous resistance is opposed more efficiently at increased $La_{m}$ , as suggested by the overshoot. This same behaviour is depicted for all the three $La$ and $Ga$ between 0.1 and 10. The system behaviour experiences a transition in the nature of its temporal response; here the transition is due to the control parameter  $La_{m}$ . In the analytical model in § 10 it is shown that this transition is indeed a standard transition between two fixed points of type node and spiral.

At a given $La_{m}$ , the increase in the $Ga$ number alters the time response in three ways. One of them has already been discussed which is that as the $Ga$ is increased, the equilibrium height for the droplet tip decreases. The second effect is that the time to reach the equilibrium state considerably reduces with increasing $Ga$ . The third observation is that at $Ga=10$ and $La_{m}=40$ , the droplet tip moves downward instead of moving upward, as expected in § 5.3. Also, provided that the value of $La_{m}$ is such that it causes an overshoot, the number of cross-overs around the steady state location increase with increasing  $Ga$ .

Whereas the droplet reaches a steady state relatively early in the case of higher $Ga$ , reaching equilibrium is apparently asymptotic in the case of $Ga=0.1$ . For $Ga=10.0$ and 1.0, the steady state is reached exactly well before $t=20$ . However when $Ga=0.1$ , the $h_{tip}$ value has not exactly reached its steady state even at $t=35$ . The droplet essentially creeps towards the equilibrium. The viscosity is thus more dominant at low $Ga$ and $La_{m}$ . This also holds true from the fact that $O(Ga\unicode[STIX]{x1D70C}^{\ast }g^{\ast })<O(\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\ast })$ as $Ga\rightarrow 0$ for $Ga<1$ and $O(La_{m}\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{m}^{\ast })\rightarrow O(\unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\ast })$ as $La_{m}\rightarrow 1$ for $La_{m}>1$ .

Figure 10. The time-dependent vertical velocity of the tip of the droplet for different $La,Ga$ and $La_{m}$ . The highest peak in $v_{tip}$ increases with $La_{m}$ .

We further elaborate on the different characteristics of the time-dependent evolution of the droplet by analysing the vertical velocity of the tip of the droplet. The results are presented in figure 10. First, the maximum rise in velocity is encountered only during the very initial transient stage. In all of the cases, starting from rest, the droplet gains a maximum velocity which again becomes negligible within a non-dimensional time of approximately $t=5.0$ . This further refines our intuition that the hydrodynamic flow is dominant during the initial rise of the droplet. It is again observable that $La$ has a negligible effect on the time-dependent evolution of the droplet and it primarily influences only the shape of the droplet. The $Ga$ and the $La_{m}$ numbers, on the other hand, considerably change the levitation velocities. The oscillations of the droplet around the final equilibrium position at increased $Ga$ are appreciably resolved via the velocity characteristics. This feature was less appreciable in the displacement plots. Now at least two periods of oscillation around the equilibrium are clearly identifiable for $Ga=10.0$ . For this value of  $Ga$ , it is depicted that the levitation velocity can also be negative after reaching its positive maximum, or even at the very beginning of the phenomenon if $La_{m}$ is low, indicating the downward motion of the droplet after the overshoot.

The temporal route to the final equilibrium state thus can be either monotonic, or through undulations. Further, the monotonic levitation can be asymptotic in nature where the droplet takes a long time to reach its final equilibrium (creep). The displacement plots have revealed the monotonicity and the asymptotic nature, while the velocity plots have resolved the undulations around the equilibrium more apparently.

Figure 11. The summary of the shapes of the droplet for different $La$ and $Ga$ at various time instants. The $La_{m}=1000$ in this case.

Figure 12. The summary for the steady state height of the tip of the droplet for different $La,Ga$ and $La_{m}$ . The steady state droplet shapes are also presented.

5.5 Summary of droplet shapes and levitation height for constant permeability

The observations under the constant ferrofluid permeability assumption are summarized in terms of the droplet shapes and the corresponding levitation height. First we discuss the shapes.

Different interfacial configurations of the droplet emerge when the combination of $La,La_{m}$ and $Ga$ is changed. These shapes also vary with time for given $La,La_{m}$ and $Ga$ . A summary of the levitating droplet shapes with changing time, $La$ and $Ga$ is presented in figure 11. The $La_{m}$ number is 1000. The droplet shapes are more stable for higher magnitude of $La$ . On the other hand, at low $La=0.1$ , highly deformed droplet shapes are observed. In this regime, the shape of the droplet is skirted near the initiation of the phenomena, at $t=0.5$ and 1.0. Afterwards, a transition to a tooth-shaped or segmented-ring-shaped configuration takes place, depending on the value of  $Ga$ . For example, for $La=0.1$ and $Ga=0.1$ , the transition is from skirted towards ring segment, while for $La=0.1$ and $Ga=10.0$ , the droplet transforms from a skirted to a tooth-like shape. These shapes are suppressed for higher $Ga$ ( $=10.0$ ). An increase in $La$ transforms the droplet shape towards an oval.

It is interesting to note that, for $La=0.1$ and 1.0, the maximum deformation of the droplet occurs between $t=1.0$ and $t=4.0$ . From the velocity plots, this is the period when the droplet is retarding after attaining a maximum velocity. However, for $La=10.0$ , the maximum deformation occurs before $t=1.0$ , when the droplet is accelerating. Thus a stiffer interface (higher  $La$ ) deforms maximally while in acceleration and a flexible interface (lower  $La$ ) exhibits maximal deformation while in retardation.

Figure 12 presents the vertical height of the droplet tip as $La_{m}$ is varied at different $La$ and $Ga$ . The levitation height is along the abscissa and the $La_{m}$ increases along the ordinate logarithmically. The three zones, separated by two dashed lines, are for different $Ga$ while in each zone, the line plots are for different  $La$ . The corresponding droplet shapes are also shown, except for the intermediate case of $La=1.0$ . A significant outcome from this semi-log graph is that the equilibrium state levitation height of the droplet increases nearly exponentially with  $La_{m}$ , for all the simulated $La$ and  $Ga$ . The roles of the non-dimensional numbers – $La,Ga$ and $La_{m}$ – are now more understandable from this plot. Whereas $Ga$ has the least effect on the droplet shape, $La$ has the least effect on the levitation height. The $La_{m}$ number influences the phenomenon significantly in both ways. The shapes in figure 12 are taken at $t=30$ . At this time the shapes are at steady state, except for the cases where the droplet creeps towards the equilibrium. However, even for the creeping cases, the shapes at $t=30$ are expected to be quite close to the equilibrium, as indicated from the displacement plots in figure 9.

7.1 The effect of parameter $\unicode[STIX]{x1D6FE}_{o}$

Apparently, the parameter $\unicode[STIX]{x1D6FE}_{o}$ has a considerable effect on the deformation of the droplet, as well as on the levitation height. For $\unicode[STIX]{x1D6FE}_{o}=0.111$ , the droplet has a nearly round shape during its rise (figure 14). A concavity starts appearing at the tail of the droplet for $\unicode[STIX]{x1D6FE}_{o}=0.278$ , and for a given time, the extent of the concavity increases as $\unicode[STIX]{x1D6FE}_{o}$ increases. The droplet shapes at $\unicode[STIX]{x1D6FE}_{o}=0.278$ and 0.389 are similar to that of a cardioid without a cusp. For $\unicode[STIX]{x1D6FE}_{o}=0.556$ , 0.833 and 1.111, the bottom feature of the droplet has widened and at $\unicode[STIX]{x1D6FE}_{o}=1.389$ and 1.667, the droplet is more flattened and skirted.

Figure 15. The effect of parameter $\unicode[STIX]{x1D6FE}_{o}$ on the field gradient across the top and bottom poles of the droplet (denoted by $\unicode[STIX]{x0394}H$ in the figure). The permeability in the expression of $La_{m}$ is now field dependent while the fixed parameters are $Ga=0.1$ , $La=0.1$ , $\unicode[STIX]{x1D709}_{o}=0.185$ . All four frames are at initial time $t=0.0$ . The field strength difference across the poles of the droplet ( $\unicode[STIX]{x0394}H$ ) increases with increasing $\unicode[STIX]{x1D6FE}_{o}$ .

In this subset of simulations, the non-dimensional numbers $La,Ga$ and $\unicode[STIX]{x1D709}_{o}$ are all kept constant. Except for the magnetic permeability, the other quantities in the expression of $La_{m}$ are also fixed. This essentially implies that the changes in the dynamics of the flow come from the non-uniform spatial distribution of  $\unicode[STIX]{x1D707}_{f}$ . This non-uniformity makes it intuitive that $\unicode[STIX]{x1D707}_{f}$ along the interface of the droplet is no longer an invariant. We seek the effects of the variability of  $\unicode[STIX]{x1D707}_{f}$ , due to different  $\unicode[STIX]{x1D6FE}_{o}$ , by plotting the contours of $H/H_{o}$ around the initial state of the droplet. The contours are depicted in figure 15. The four frames are for $\unicode[STIX]{x1D6FE}_{o}=1.667$ , 1.111, 0.556 and 0.111. The spatial non-uniformity of $H/H_{o}$ is sensitive to the parameter  $\unicode[STIX]{x1D6FE}_{o}$ . The difference in the field strengths at the top and the bottom pole of the droplet ( $\unicode[STIX]{x0394}H$ ), or in other words, the field gradient along the vertical axis of the droplet, in a global sense, increases with increasing  $\unicode[STIX]{x1D6FE}_{o}$ . The local gradients across this axis are relatively smaller for lower values of  $\unicode[STIX]{x1D6FE}_{o}$ , for example in the plot for $\unicode[STIX]{x1D6FE}_{o}=0.111$ . Thus, a rational explanation is that the increase in $\unicode[STIX]{x1D6FE}_{o}$ increases the field gradient across the droplet, and thus the acceleration of the levitating droplet is higher for larger $\unicode[STIX]{x1D6FE}_{o}$ (figures 14 and 15).

The fundamental change in the droplet response due to the nonlinearity caused by the additional parameter $\unicode[STIX]{x1D6FE}_{o}=(3\unicode[STIX]{x1D712}_{o}\unicode[STIX]{x1D709}_{o})$ , which basically controls the magnetization curve of the ferrofluid, is deducible by deriving asymptotic expressions for the order of the magnetic force in the limits $\unicode[STIX]{x1D6FE}_{o}\rightarrow 0$ and $\unicode[STIX]{x1D6FE}_{o}\rightarrow \infty$ . By transforming the following expression for the field-dependent permeability

(7.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{f}=\unicode[STIX]{x1D707}_{o}(1+\unicode[STIX]{x1D712}(H))=\unicode[STIX]{x1D707}_{o}\left(1+\frac{M_{s}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\left[\coth \unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|-\frac{1}{\unicode[STIX]{x1D6FE}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\right]\right),\end{eqnarray}$$

to its non-dimensional form

(7.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{f}=\unicode[STIX]{x1D707}_{o}\left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|}\left[\coth \unicode[STIX]{x1D6FE}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|-\frac{1}{\unicode[STIX]{x1D6FE}_{o}|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|}\right]\right),\end{eqnarray}$$

that for bounded $|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}^{\ast }|$ ,

(7.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D707}_{f}\sim \unicode[STIX]{x1D707}_{o},\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow 0,\quad \text{and,}\\ \displaystyle \unicode[STIX]{x1D707}_{f}\sim \unicode[STIX]{x1D707}_{o}\left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}}\right),\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow \infty .\end{array}\right\}\end{eqnarray}$$

Using the above asymptotic relations, the order of the magnetic force, say  $f_{m}$ , becomes

(7.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle f_{m}\sim \left(\frac{\unicode[STIX]{x1D707}_{f}H_{o}^{2}}{R}\right)\sim \left(\frac{\unicode[STIX]{x1D707}_{o}H_{o}^{2}}{R}\right),\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow 0,\quad \text{and,}\\ \displaystyle f_{m}\sim \left(\frac{\unicode[STIX]{x1D707}_{f}H_{o}^{2}}{R}\right)\sim \left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}}\right)\left(\frac{\unicode[STIX]{x1D707}_{o}H_{o}^{2}}{R}\right),\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow \infty .\end{array}\right\}\end{eqnarray}$$

Rewriting the quantity $\unicode[STIX]{x1D707}_{o}H_{o}^{2}/R$ as $\unicode[STIX]{x1D709}_{o}^{2}\unicode[STIX]{x1D707}_{o}M_{s}^{2}/R$ (as $H_{o}=\unicode[STIX]{x1D709}_{o}M_{s}$ ), the order of the magnetic force is rewritten to

(7.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle f_{m}\sim \unicode[STIX]{x1D709}_{o}^{2}\left(\frac{\unicode[STIX]{x1D707}_{o}M_{s}^{2}}{R}\right),\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow 0,\quad \text{and,}\\ \displaystyle f_{m}\sim \unicode[STIX]{x1D709}_{o}^{2}\left(1+\frac{1}{\unicode[STIX]{x1D709}_{o}}\right)\left(\frac{\unicode[STIX]{x1D707}_{o}M_{s}^{2}}{R}\right),\quad \text{as}~\unicode[STIX]{x1D6FE}_{o}\rightarrow \infty .\end{array}\right\}\end{eqnarray}$$

The factor in the expression for the order of magnetic force is $\unicode[STIX]{x1D709}_{o}^{2}(1+1/\unicode[STIX]{x1D709}_{o})$ . The value of this function has a global minimum at $\unicode[STIX]{x1D709}_{o}=-1/2$ . As by definition $H_{o}$ and $M_{s}$ are positive, $\unicode[STIX]{x1D709}_{o}$ is also positive in the interval of interest, and thus any increase in $\unicode[STIX]{x1D709}_{o}$ (or equivalently  $\unicode[STIX]{x1D6FE}_{o}$ ) will also increase  $f_{m}$ . However, as the behaviour of $f_{m}$ near $\unicode[STIX]{x1D6FE}_{o}=0$ and in the limit $\unicode[STIX]{x1D6FE}_{o}\rightarrow \infty$ is different, its effect on the droplet response in the two regimes is also expected to be different. This is depicted by simulations in figure 14 where the droplet responds differently with increasing  $\unicode[STIX]{x1D6FE}_{o}$ .

8.1 Single and multiple stable and unstable states

First, we study the field contours for three different arrangements of the magnetic poles at the boundaries. The three arrangements are (i) all of the four walls are covered with alternating magnetic poles, (ii) the bottom and side walls are covered by alternating poles of magnets and only the top wall is free to the field and (iii) only the bottom wall is covered by the magnetic poles while the top and the side walls are free. Figure 17 shows the absolute field contours around the initial state of the eccentrically placed droplet for these three cases.

Figure 17. The absolute magnetic field contours around the eccentrically placed droplet when – from left to right – (a) all of the four walls are covered with alternate magnetic poles, (b) the bottom and both the side walls are covered and (c) only the bottom wall is covered.

When all of the four walls are covered with alternating magnetic poles, a minimum field region is created near the centre of the domain (figure 17 a). The field gradient is everywhere pointing inward to the domain and it is expected to be the equilibrium location for the droplet after levitation.

In the second, particularly interesting case (figure 17 b), where only the top wall is free from the magnetic poles, two minimum field regions are established, which indicates the possibility of multiple stable states. The two low field regions give rise to a cat-eye-like field configuration. As the droplet is closer to the right eye, there is an increased probability that it will be trapped in that region.

In figure 17(c), the top and the side walls are now free from the magnetic poles while the bottom wall field boundary condition is the same as before. In this case, the droplet is expected to levitate, but as the field gradient is diverging or equivalently the contours are concave downward throughout the domain, the stability may not be ensured.

Figure 18. The interface of the eccentrically placed droplet when – from left to right – (a) all of the four walls are covered with alternate magnetic poles, (b) the bottom and both the side walls are covered and (c) only the bottom wall is covered.

The movement of the droplet is now simulated and the validity of the above intuitions, made by looking at the initial field contours alone, is tested. The movement of the droplet under the above stated three magnet arrangements is shown in figure 18. In the first case the droplet seeks a stable horizontal location, $x=0.5$ . The vertical location is such that the magnetic forces balance the weight of the droplet. In the second case, although there are two minimum field locations, the droplet seeks the nearby one, as argued previously. In the third and the final case, the path of the droplet levitation is unstable. The equilibrium in this case is only possible when a physical support from the side walls is present; otherwise, it can only be marginally stable if a perfect symmetry is maintained during the experiment. The mechanism is further explored below.

8.2 Interplay between magnetic forces and vorticity: stability of levitation path

One non-trivial and counter-intuitive observation is made in regard to the unstable path of the droplet levitation in figure 18(c). At the first look, the field contours in figure 17(c) are such that their normal is pointing in the direction shown by the dashed line. As this direction is parallel to $\unicode[STIX]{x1D735}H$ , the droplet is expected to move along this direction. However, in figure 18(c), the path of the droplet exhibits an opposite behaviour. For some initial time, the droplet moves along $\unicode[STIX]{x1D735}H$ , alters its direction and then starts moving towards the north-west. This behaviour of the droplet levitation directs that although the instability of the levitation is predictable from the configuration of the absolute magnetic field, the predictability of the actual unstable path is more involved. Essentially, the information about the magnetic field patterns alone is not sufficient to predict the path of the droplet if the levitation becomes unstable. The magnetic field contours and the hydrodynamic flow field around the droplet should be analysed in a coupled way.

Figure 19. The maximum value of the normalized absolute vorticity (a), and the velocity vector field (b) in and around the levitating droplet. The initial horizontal location of the droplet is central in this case.

We study the velocity field solution when the eccentricity in the horizontal location of the droplet is zero. The same is presented in figure 19. In this case, there is no sideways movement of the droplet and the path of levitation is a straight vertical line. The magnetic field is due to only two magnets with alternating poles at the bottom wall, while the other three walls are not covered with the magnets. As can be observed, two counter-rotating vortices originate at the tail of the droplet during the very initiation of the phenomenon. The situation is horizontally symmetric and there are no net lateral forces on the droplet due to the flow dynamics. The flow, however, does not show a similar character when the initial location of the droplet is not horizontally central. The velocity field solution for this situation is shown in figure 20. The asymmetry in this situation causes one of the vortices, which is closer to the horizontal centre ( $x=0.5$ ) to have more size and strength. This is evident by looking the left counter-clockwise rotating vortex in figure 20 for $t=0.2$ , and the corresponding vorticity map. The magnetic forces have given rise to a vortex structure which is actually tilting the droplet laterally away from the expected path of the droplet. The levitation path eventually deviates from the direction along $\unicode[STIX]{x1D735}H$ due to a coupled interplay between the magnetic and the hydrodynamic flow field.

Figure 20. The maximum value of the normalized absolute vorticity (a), and the velocity vector field (b) in and around the levitating droplet. The initial horizontal location of the droplet is eccentric in this case.

The arrangement of the magnets has a vital role in the stability of the levitation. In the absence of the magnets at the top and side walls, the lateral perturbations in the droplet path do not decay. The selection of the magnetic source arrangement, thus, should be a carefully considered design feature and stable levitation is certainly possible under an appropriately chosen magnetic field.

Based on the above observations from the simulations, a schematic can be proposed in terms of the initial field contour configurations, which can help to depict the stability, instability or multiple stabilities of the droplet levitation. The schematic is presented in figure 21. The mechanism is analogous to the classical ball-on-ramp example which is often used to understand the concept of stability, and this analogy can be applied to study the stability of levitation by solving for the absolute field contours in the domain produced by a given arrangement of magnets. The stability of the droplet path, however, can be non-trivial and should be analysed in conjunction with the hydrodynamic flow field near the droplet.

Figure 21. The schematic for determining the stability of the levitation: (a) stable, (b) unstable and (c) possibility of multiple stable states. $\unicode[STIX]{x0394}H$ is the difference in the field at the bottom and the top of the droplet.

9.1 Singularity at the tail of the droplet

The shapes of the droplet at different instants of time from the experiment and simulation are compared in figure 24. For the above characterized samples of fluids, we notice that the levitation height is small, which is due to the weakly magnetizable nature of the ferrofluid sample ( $M_{s}=57.7,\unicode[STIX]{x1D712}_{o}=0.0189$ ). Interestingly, a formation of singular projection at the tail of the levitating droplet is observed in both the experiment and the simulation. The central portion of the tail of the droplet remains adhered in the vicinity of the bottom wall for some initial time, although there is a layer of ferrofluid separating it from the bottom wall. The width of this adhered part of the interface reduces gradually and, at approximately $t=0.035~\text{s}$ in the experiment, it detaches from the bottom wall. The singularity is formed after this instant at the tail of the droplet ( $t=0.036~\text{s}$ ). As the droplet levitates, the singular point appears to being pulled towards the bulk of the droplet. Eventually the singularity disappears and the droplet tail becomes smooth and concave downward. The simulations have predicted this phenomenon reasonably well. This also strongly supports our argument in § 6 that the interfacial singularities or cusps at the interface of the levitating droplet in a ferrofluid can be ubiquitous, especially for low interfacial tension (low  $La$ ).

Although the capture of the singular shape requires full simulations, a scaling argument for the interface shape just before the appearance of the singularity is given as follows. As initially a segment of the droplet interface is settled near the bottom wall, the field strength along this portion of the interface is approximated as $H(x)|_{y=0}=-H_{o}\sin (2\unicode[STIX]{x03C0}x/L)$ , according to the field boundary condition at the bottom wall. For a small initial time $\unicode[STIX]{x1D6FF}t$ we assume that the magnetic field at this segment of the interface is close to the field value at the wall itself. The order of the vertical magnetic force on this segment of the interface is then $f_{m}\sim \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{m}\sim -(1/2)H^{2}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ , where $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ is the permeability difference across the interface. Using the field strength along this interface segment, $f_{m}\sim -(1/2)H_{o}^{2}\sin ^{2}(2\unicode[STIX]{x03C0}x/L)\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ . For a small initial period of time $\unicode[STIX]{x1D6FF}t$ , the segment will rise faster if $f_{m}$ at the wall is higher and vice versa. If the displaced vertical location of this segment after $\unicode[STIX]{x1D6FF}t$ is $\unicode[STIX]{x1D6FF}y$ , then $\unicode[STIX]{x1D6FF}y(x)\sim f_{m}(x)\sim \sin ^{2}(x)$ . The squared sine interface description at the tail of the droplet is closely resembled at $t=0.030~\text{s}$ in the simulations and at $t=0.035~\text{s}$ in the experiments in figure 24, where the tip of the projection is still close to the wall and is not cusped (also indicated by smoothness of $\sin ^{2}(x)$ ). However, at relatively larger time the $\sin ^{2}(x)$ feature transitions towards a cusp, and the approximation that the field value at the segment is equal to the field value at the bottom wall is no more valid. After this regime, the hydrodynamics and the effect of higher outer liquid viscosity (see § 6) become significant.

The sharp singular tips are observed in droplet systems, for example in Stone, Lister & Brenner (Reference Stone, Lister and Brenner1999) under an electric field and by Rowghanian et al. (Reference Rowghanian, Meinhart and Campàs2016) for a magnetic drop under a uniform magnetic field. Stone et al. (Reference Stone, Lister and Brenner1999) have concluded that the conical tips appear for a certain threshold electric field strength and dielectric constant ratio. Rowghanian et al. (Reference Rowghanian, Meinhart and Campàs2016) focus on curved tip regimes, and clearly discuss that sharp tips are possible due to destabilization of the droplet under high field strengths and high localization of the magnetic stresses at the droplet poles. In the present study, the situation is, however, dynamic and the field is non-uniform. A low interfacial tension (low Laplace number  $La$ ) in combination with the hydrodynamic viscous stress difference (exhibited through viscosity ratio change) across the interface are two additional important factors which can help realize such sharp conical tips. In practical situations, the former is reasonably modifiable using a surfactant solution and the latter by taking a different non-magnetic droplet medium in a ferrofluid (or a non-magnetic medium outside of a ferrofluid droplet).

Although the droplet is essentially moving in a quasi-two-dimensional cell, and it might be more appropriate to use gap-averaged Navier–Stokes equations, it turns out that the current two-dimensional mathematical description captures the droplet’s shape reasonably well. Thus the spatial solution is quite accurate. However, during the later stages of the droplet rise, the solution is temporally deviating (figure 24). As far as the spatial shapes of the droplet are concerned, any three-dimensional effects due to the out-of-plane curvature of the droplet in the Hele-Shaw gap are not dominant. We attribute the accurate capturing of the droplet interface shape to the following points: (i) careful characterization of the ferrofluid sample nonlinear magnetization curve, as well as some other properties, and implementing of the same in the simulations using Langevin’s nonlinear magnetization model (figure 22), (ii) implementing the initial interface condition precisely similar to that which we find for a resting droplet in the experiments, (iii) convenience in closely realizing the actual field boundary conditions due to the Halbach array, (iv) modelling similar flow boundary conditions as in the experiments and (v) using state-of-the-art numerical techniques for computational accuracy. In addition, the ferrofluid sample in the experiments is characterized and found to be weakly magnetizable. The numerical solution is then expected to be more accurate due to a less sharp change in the magnetic properties and field (and thus a softer jump in the magnetic force term in the Navier–Stokes equations) across the interface. The possible reasons why temporally the two-dimensional assumption is only reasonable might be due to the following unrelaxed assumptions: (i) neglecting the out-of-plane curvature in simulations due to the Hele-Shaw geometry, or in other words, neglecting the quasi-two-dimensional nature of the problem and considering it as a two-dimensional problem in the simulations, (ii) neglecting an important aspect of the magnetization relaxation dynamics in ferrofluids and considering instantly relaxing magnetization when the magnetic field changes. We believe that incorporation of the magnetization relaxation time effect coupled with a two-phase model for the ferrofluid interface, together with corrections due to the third dimension (e.g. implementing the Darcy law (Gondret & Rabaud Reference Gondret and Rabaud1997) together with inertial corrections for it (Ruyer-Quil Reference Ruyer-Quil2001)), can closely resolve the temporal details.

10.1 Levitation height and onset condition

The differential equation (10.4) is a second-order ordinary but nonlinear equation due to the form of $\unicode[STIX]{x1D6F9}(\tilde{\unicode[STIX]{x1D701}})$ . It is reduced to a first-order system of two differential equations by introducing a variable $\unicode[STIX]{x1D6F1}=\dot{\tilde{\unicode[STIX]{x1D701}}}$ , written as

(10.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\dot{\tilde{\unicode[STIX]{x1D701}}}=\unicode[STIX]{x1D6F1},\\ \dot{\unicode[STIX]{x1D6F1}}=-{\mathcal{V}}\unicode[STIX]{x1D6F1}-{\mathcal{G}}+\unicode[STIX]{x1D6F9}(\tilde{\unicode[STIX]{x1D701}}).\end{array}\right\}\end{eqnarray}$$

The long-time equilibrium or fixed point of the system is obtained by setting the droplet speed $\dot{\tilde{\unicode[STIX]{x1D701}}}=\unicode[STIX]{x1D6F1}$ , and acceleration $\dot{\unicode[STIX]{x1D6F1}}$ to 0 in (10.5), which for non-zero  $\unicode[STIX]{x1D6F9}$ , gives the steady state fixed point $(\unicode[STIX]{x1D6F1}_{\ast },\tilde{\unicode[STIX]{x1D701}}_{\ast })$

(10.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6F1}_{\ast }=0,\quad \text{steady state velocity},\\ \displaystyle \unicode[STIX]{x1D6F9}(\tilde{\unicode[STIX]{x1D701}}_{\ast })={\mathcal{G}},\quad \text{equation for steady state levitation height}.\end{array}\right\}\end{eqnarray}$$

As $\unicode[STIX]{x1D6F9}(\tilde{\unicode[STIX]{x1D701}})={\mathcal{L}}\text{e}^{-{\mathcal{K}}\tilde{\unicode[STIX]{x1D701}}}$ , it turns out that

(10.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6F1}_{\ast }=0,\quad \text{steady state velocity},\\ \displaystyle \tilde{\unicode[STIX]{x1D701}}_{\ast }=-\frac{1}{{\mathcal{K}}}\ln \frac{{\mathcal{G}}}{{\mathcal{L}}}=\frac{1}{{\mathcal{K}}}\ln \frac{{\mathcal{L}}}{{\mathcal{G}}},\quad \text{steady state levitation height}.\end{array}\right\}\end{eqnarray}$$

The steady state levitation height of the droplet depends on parameter ${\mathcal{K}}$ and the ratio ${\mathcal{G}}/{\mathcal{L}}$ (which is $\propto Ga/La_{m}$ ), and is independent of the parameter  ${\mathcal{V}}$ . This is expected as ${\mathcal{V}}$ is the strength of the viscous drag force opposing the droplet motion and it must vanish at steady state. If the levitation height $\tilde{\unicode[STIX]{x1D701}}_{\ast }$ is to be a non-zero positive, then $(1/{\mathcal{K}})\ln ({\mathcal{L}}/{\mathcal{G}})\geqslant 0$ , which provides the condition for the onset of levitation

(10.8a,b ) $$\begin{eqnarray}{\mathcal{L}}\geqslant {\mathcal{G}},\quad \text{or},\quad \frac{La_{m}}{Ga}\geqslant \unicode[STIX]{x1D6FC}_{1}\frac{(1-\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D70C}_{d})}{(1-\unicode[STIX]{x1D707}_{o}/\unicode[STIX]{x1D707}_{f})},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{1}=1/{\mathcal{K}}$ is constant. If the onset condition is satisfied, the levitation height is

(10.9) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D701}}_{\ast }=\unicode[STIX]{x1D6FC}_{1}\ln \frac{La_{m}(1-\unicode[STIX]{x1D707}_{o}/\unicode[STIX]{x1D707}_{f})}{Ga(1-\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D70C}_{d})}.\end{eqnarray}$$

The above expression is physically is in agreement with the facts – the rise of the droplet will increase in a ferrofluid sample of higher permeability, with lower droplet mass densities, with the increase in $La_{m}$ or decrease in  $Ga$ .

10.2 Transition in the nature of the fixed point

The transitions in the nature of the fixed point are revealed by evaluating the eigenvalues of the Jacobian $\unicode[STIX]{x1D645}$ for the system (10.5), which is

(10.10) $$\begin{eqnarray}\unicode[STIX]{x1D645}=\left[\begin{array}{@{}cc@{}}0 & 1\\ \displaystyle \frac{\text{d}\unicode[STIX]{x1D6F9}}{\text{d}\tilde{\unicode[STIX]{x1D701}}} & -{\mathcal{V}}\end{array}\right].\end{eqnarray}$$

If $\unicode[STIX]{x1D706}$ are the eigenvalues of $\unicode[STIX]{x1D645}$ , then the characteristic equation $|\unicode[STIX]{x1D645}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D644}|=0$ is $\unicode[STIX]{x1D706}^{2}+{\mathcal{V}}\unicode[STIX]{x1D706}-\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}}=0$ and the eigenvalues or the characteristic values of the fixed point, $\unicode[STIX]{x1D706}_{1,2}$ , are $\unicode[STIX]{x1D706}_{1,2}=-({\mathcal{V}}/2)\pm (1/2)({\mathcal{V}}^{2}-4(-\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}}))^{1/2}$ . This gives the possibilities ${\mathcal{V}}^{2}>-4(\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}})$ (distinct real eigenvalues), ${\mathcal{V}}^{2}=-4(\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}})$ (identical real eigenvalues) and ${\mathcal{V}}^{2}<-4(\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}})$ (complex conjugate eigenvalues). The first case physically suggests that the fixed point is of type node and the droplet will approach the steady state monotonically. The third case of complex eigenvalues, however, suggest that the fixed point is of type spiral and the droplet will oscillate around the fixed point before reaching the steady state. The second case gives the critical value of the parameter  ${\mathcal{V}}$ , at which the stable fixed point $\unicode[STIX]{x1D701}_{\ast }$ changes its type from spiral to node or vice versa, according to

(10.11) $$\begin{eqnarray}\displaystyle {\mathcal{V}}_{c}^{2}=-4\left[\frac{\text{d}\unicode[STIX]{x1D6F9}}{\text{d}\tilde{\unicode[STIX]{x1D701}}}\right]_{\tilde{\unicode[STIX]{x1D701}}_{\ast }}=-4\left[\frac{\text{d}}{\text{d}\tilde{\unicode[STIX]{x1D701}}}{\mathcal{L}}\text{e}^{-{\mathcal{K}}\tilde{\unicode[STIX]{x1D701}}}\right]_{\tilde{\unicode[STIX]{x1D701}}_{\ast }}=4{\mathcal{K}}{\mathcal{G}}.\end{eqnarray}$$

Thus from the condition of complex conjugate eigenvalues, ${\mathcal{V}}^{2}<-4(\text{d}\unicode[STIX]{x1D6F9}/\text{d}\tilde{\unicode[STIX]{x1D701}})$ , oscillations around the equilibrium point do occur if

(10.12a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{V}}^{2}<4{\mathcal{K}}{\mathcal{G}},\quad \text{or},\quad \frac{\unicode[STIX]{x1D702}_{f}}{\unicode[STIX]{x1D702}_{d}}<\unicode[STIX]{x1D6FC}_{2}\sqrt{\left(1-\frac{\unicode[STIX]{x1D70C}_{f}}{\unicode[STIX]{x1D70C}_{d}}\right)Ga},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{2}=(4/9){\mathcal{K}}^{1/2}$ is constant. It is important to note that the parameter ${\mathcal{L}}$ does not influence the condition for the oscillations around the equilibrium; it only appears in the expression of the steady state levitation height.

Figure 25. The behaviour of the solution of (10.3) near the equilibrium/fixed point, depicted by the phase portraits in the displacement–velocity $(\tilde{\unicode[STIX]{x1D701}},\dot{\tilde{\unicode[STIX]{x1D701}}})$ plane with varying ${\mathcal{V}}$ ( ${\mathcal{V}}$ is proportional to the viscosity ratio). The nature of the fixed point $(\tilde{\unicode[STIX]{x1D701}}_{\ast },\dot{\tilde{\unicode[STIX]{x1D701}}}_{\ast })=((1/{\mathcal{K}})\ln ({\mathcal{L}}/{\mathcal{G}}),0)$ changes at critical ${\mathcal{V}}_{c}=1.804$ from a stable spiral to a stable node. Here, ${\mathcal{L}}/{\mathcal{G}}=1.21$ and ${\mathcal{K}}=\unicode[STIX]{x03C0}/5$ .

Figure 26. The comparison of the time displacement of the droplet between the simulations and the analytical model for two extreme cases of $Ga$ ( $=0.1$ and 10). The $La_{m}$ for the two cases is chosen such that the levitation height $\tilde{\unicode[STIX]{x1D701}}_{\ast }$ ( $\propto La_{m}/Ga$ ) remains in the same range. Both approaches have predicted that there is a transition between a node and a spiral.

The nature of the stable fixed point $\unicode[STIX]{x1D701}_{\ast }$ and transitions in its behaviour are further studied through phase portraits in the $(\unicode[STIX]{x1D701},\dot{\unicode[STIX]{x1D701}})$ plane. The differential equation for $(\unicode[STIX]{x1D701},\dot{\unicode[STIX]{x1D701}})$ from (10.5) is

(10.13) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}\tilde{\unicode[STIX]{x1D701}}}=\frac{\text{d}\dot{\tilde{\unicode[STIX]{x1D701}}}}{\text{d}\tilde{\unicode[STIX]{x1D701}}}=\frac{\dot{\unicode[STIX]{x1D6F1}}}{\dot{\tilde{\unicode[STIX]{x1D701}}}}=\frac{-{\mathcal{V}}\unicode[STIX]{x1D6F1}-{\mathcal{G}}+\unicode[STIX]{x1D6F9}(\tilde{\unicode[STIX]{x1D701}})}{\unicode[STIX]{x1D6F1}}.\end{eqnarray}$$

Using this, first we show the effect of variation of ${\mathcal{V}}$ ( $\propto$ viscosity ratio) at fixed ${\mathcal{G}}$ ( $\propto Ga$ ) and ${\mathcal{L}}$ ( $\propto La_{m}$ ). This is depicted in figure 25 for different initial conditions $(\unicode[STIX]{x1D701}(0),\dot{\unicode[STIX]{x1D701}}(0))$ on a circle around the fixed point. For $\unicode[STIX]{x1D702}_{f}/\unicode[STIX]{x1D702}_{d}=0.5$ , ${\mathcal{V}}=2.25$ , and thus we vary ${\mathcal{V}}$ around $O(1)$ . The behaviour of the linearized solution around the fixed point $\unicode[STIX]{x1D701}_{\ast }$ is clearly sensitive to  ${\mathcal{V}}$ , although the actual magnitude of $\unicode[STIX]{x1D701}_{\ast }$ depends on ${\mathcal{L}}/{\mathcal{G}}$ . The variation in ${\mathcal{V}}$ changes the nature of the fixed point from being a spiral to a monotonic attractor at  ${\mathcal{V}}_{c}$ . This observation is in good agreement with the simulations in figures 9 and 10, where the transition between monotonicity and oscillations of the curves was noted.

We conclude by explicitly comparing the solution of the analytical model (10.3) to the simulations of figure 9. In the simulations we noted that the response changes from monotonic to undulatory when $Ga$ is increased (figure 9). Here we take the extreme cases of $Ga=0.1$ and 10 considered there. The comparison is shown in figure 26. The corresponding trajectories in the ( $\tilde{\unicode[STIX]{x1D701}},\dot{\tilde{\unicode[STIX]{x1D701}}}$ ) plane are also depicted in the inset. The $La_{m}$ value is chosen from the results of figure 9 such that the ratio ${\mathcal{L}}/{\mathcal{G}}$ is more or less the same in the dynamical model. This results in nearly the same final levitation height according to $\tilde{\unicode[STIX]{x1D701}}_{\ast }={\mathcal{K}}^{-1}\ln ({\mathcal{L}}/{\mathcal{G}})$ , and thus only the nature of the response can be closely compared. We have obtained a reasonable agreement between the simulations and the dynamical model regarding the statement of transition between monotonic and undulatory droplet levitation.

Thus the model captures the essential features of the levitation dynamics of the droplet, and at the same time, it provides a simpler alternative in place of fully fledged simulations to predict the condition for the onset of levitation, the levitation height, the condition for oscillations around the equilibrium and the nature of the solution around the equilibrium. The model, however, does not care about the lateral stability of the levitation path. Those conditions for the lateral stability need a coupled analysis of the magnetic field and the flow solutions (§ 8).