Tuning a magnetic field to generate spinning ferrofluid droplets with controllable speed via nonlinear periodic interfacial waves
TTuning a Magnetic Field to Generate Spinning Ferrofluid Droplets with ControllableSpeed via Nonlinear Periodic Interfacial Waves
Zongxin Yu ∗ and Ivan C. Christov † School of Mechanical Engineering, Purdue University, West Lafayette, Indiana 47907, USA (Dated: September 11, 2020)Two dimensional free surface flows in Hele-Shaw configurations are a fertile ground for exploringnonlinear physics. Since Saffman and Taylor’s work on linear instability of fluid–fluid interfaces,significant effort has been expended to determining the physics and forcing that set the linear growthrate. However, linear stability does not always imply nonlinear stability. We demonstrate how thecombination of a radial and azimuthal external magnetic field can manipulate the interfacial shapeof a linearly unstable ferrofluid droplet in a Hele-Shaw configuration. We show that weakly nonlineartheory can be used to tune the initial unstable growth. Then, nonlinearity arrests the instability,and leads to a permanent deformed droplet shape. Specifically, we show that the deformed dropletcan be set into motion with a predictable rotation speed, demonstrating nonlinear traveling waveson a fluid-fluid interface. The most linearly unstable wavenumber and the combined strength of theapplied external magnetic fields determine the traveling wave shape, which can be asymmetric.
Introduction
Recently, there has been significant in-terest in the physics of active and responsive fluids [1, 2].For example, swimming bacteria can take a suspensionof microscopic gears out of equilibrium and extract rec-tified (useful) work out of an otherwise random system[3]. One promising approach to creating active fluidswith controllable properties and behaviors is by suspend-ing many mechanical microswimmers made from shape-programmable materials [4] and actuating them with anexternal magnetic field [5, 6]. This actuation mecha-nism is particularly enticing for biological applicationsdue to the safe operation of magnetic fields in the med-ical setting (for, e.g. , targeted therapies and drug de-livery in vivo ) [7]. Even simpler than a suspension ofmagnetically-responsive mechanical microswimmers is asuspension of ferrofluid droplets, which can also respondto external magnetic fields [8, 9]. However, the ferrofluiddroplet’s interface motion and response to different typesof external magnetic fields is not well understood. Previ-ous work has addressed the linear stability of such fluid–fluid interfaces [10, 11], including stationary shapes [12],but not a droplet’s nonlinear dynamics or controllablemotion. Guided by the well-established ability of nonlin-earity to “arrest” long-wave instabilities [13], we demon-strate, using theory and nonlinear simulation, that it ispossible to “grow” linearly unstable ferrofluid interfacesinto well-defined permanent shapes. The resulting coher-ent droplet shapes are predictable and controllable via anexternal magnetic field. Importantly, these droplets canbe set into rotational motion with velocities predictableby the proposed theory, leading to the possibility of anexternally-actuated active fluid suspension.
Governing equations
We study the dynamics of aninitially circular ferrofluid droplet (radius R ) confined ina Hele-Shaw cell with gap thickness b and surroundedby air (negligible viscosity), as shown in Fig. 1, because“[i]f any [ferro]fluid mechanics problem is likely to beaccessible to theory and to direct comparison of theory and experiment it should be this one” [14]. Both flu-ids are considered incompressible. We propose to ap-ply the radially-varying external magnetic field H = H L r ˆe r + I πr ˆe θ . Anti-Helmholtz coils produce the radialcomponent, where H is a constant and L is a lengthscale [12, 15]. A long wire through the origin, carryingan electric current I , produces the azimuthal component.The combined magnetic field H forms an angle with theinitially undisturbed interface [16]. The droplet experi-ences a body force ∝ | M | ∇ | H | , where M is the mag-netization. To study shape dynamics, we assume theferrofluid is uniformly magnetized, M = χ H , where χ isits constant magnetic susceptibility [17].Enforcing no-slip on the confining boundaries and ne-glecting inertial terms, the confined flow is governed bya modified Darcy’s law [15] with gap-averaged velocity: v = − b η ∇ ( p − Ψ) , ∇ · v = 0 , (1)where p is the pressure in the droplet, η is the ferrofluid’sviscosity, Ψ = µ χ | H | / µ is the free-space per-meability. At the boundary of the droplet, the pressureis given by a modified Young–Laplace law [8, 9]: p = τ κ − µ M · ˆn ) , (2)where τ is the constant surface tension, and κ is the cur-vature of the droplet shape. The second term on theright-hand side of Eq. (2) is the magnetic normal trac-tion [8, 9], where ˆn denotes the outward unit normal vec-tor at the interface. This contribution breaks the sym-metry of the initial droplet interface, due to the projec-tion of H onto ˆn , and causes the droplet to rotate. Thekinematic boundary condition v n = − b η ∇ ( p − Ψ) · ˆn requires that the droplet boundary is a material surface. Mathematical analysis
We employ the weakly nonlin-ear approach [18] previously adapted to ferrofluid inter- a r X i v : . [ n li n . PS ] S e p FIG. 1. Schematic illustration of a Hele-Shaw cell confin-ing a ferrofluid droplet, initially circular with radius R . Theazimuthal magnetic field H a is produced by a long wire con-veying an electric current I . The radial magnetic field H r is produced by a pair of anti-Helmholtz coils with equal cur-rents I AH in opposite directions. The combined external mag-netic field H deforms the droplet, and its interface is given by h ( θ, t ). In comparison, the fluid outer to the droplet ( e.g. ,air) is assumed to have negligible viscosity and velocity. facial dynamics ( e.g. , [12, 15, 16]). The droplet inter-face is written as h ( θ, t ) = R + ξ ( θ, t ), where ξ ( θ, t ) = (cid:80) + ∞ k = −∞ ξ k ( t ) e ikθ represents the perturbation of the ini-tially circular interface, with complex Fourier amplitudes ξ k ( t ) ∈ C and azimuthal wavenumbers k ∈ Z . The veloc-ity potential φ = p − Ψ is then expanded into a Fourierseries φ ( r, θ, t ) = (cid:80) k (cid:54) =0 φ k ( t ) (cid:0) rR (cid:1) | k | e ikθ , and φ k is ex-pressed in terms of ξ k through the kinematic boundarycondition. Substituting these relations and Eq. (2) intoEq. (1), keeping only terms up to second order in ξ , wefind the dimensionless equations of motion ( k (cid:54) = 0) [19]:˙ ξ k = Λ( k ) ξ k + (cid:88) k (cid:48) (cid:54) =0 F ( k, k (cid:48) ) ξ k (cid:48) ξ k − k (cid:48) + G ( k, k (cid:48) ) ˙ ξ k (cid:48) ξ k − k (cid:48) . (3)From mass conservation, ξ = − (cid:80) k> | ξ k | /R ∀ t ≥ k ) = | k | R (1 − k ) (cid:124) (cid:123)(cid:122) (cid:125) surface tension − Ba R | k | + 2(1 + χ )N Br | k |− χ √ N Ba N Br R ik | k | (4)denotes the (complex) linear growth rate, andN Ba = µ χI π τ L , N Br = µ χH L τ (5)are the magnetic Bond numbers quantifying the ratio ofazimuthal and radial magnetic forces to capillary forces,respectively. Terms multiplied by χ arise from the normalmagnetic stress. The time and length scales used are12 ηL /τ b and L , respectively. Linear regime
First, consider Eq. (3), neglectingquadratic terms in ξ k , then (cid:60) [Λ( k )] = λ ( k ) governs theexponential growth or decay of infinitesimal perturba-tions. For λ ( k ) >
0, the interface is unstable. Specif-ically, Eq. (4) indicates that the radial magnetic fieldterm ∝ (1 + χ )N Br is destabilizing, while the azimuthalterm ∝ N Ba and surface tension are stabilizing. The mostunstable mode k m solves dλ ( k ) /dk = 0: k m = (cid:115) (cid:20) − Ba R + 2(1 + χ )N Br R (cid:21) . (6)It characterizes the dominant (cid:98) k m (cid:99) -fold symmetry of apattern. Note that the normal stress from the azimuthalmagnetic field does not contribute to the linear dynamics.The phase velocity of each mode, v p = −(cid:61) [Λ( k )] /k =2 χ √ N Ba N Br k/R in the linear regime, is set by (cid:61) [Λ( k )].A periodic shape on [0 , π ] forms a closed curve, meaningwave propagation is manifested as rotation of the droplet.Motion is caused by the magnetic normal stresses arisingfrom the combined magnetic field [20]. This linear anal-ysis indicates that perturbations of the droplet interfacecan propagate (and, since v p = v p ( k ), they also expe-rience dispersion). Such wavepackets will either decayor blow-up exponentially according to the sign of λ ( k ).However, this is not the whole story, and nonlinearly sta-ble traveling shapes exist, as we now show. Nonlinear regime
To demonstrate the possibility ofnonlinear traveling waves in this system, we numericallysolve the weakly nonlinear mode-coupling equations (3)for five modes ( i.e. , k, k, . . . , k ). The fundamentalmode k = 7 is chosen to allow propagating solutions overa wider swath of the (N Ba , N Br , k m ) space (compared tochoosing k < k > | c n | and phases ∠ [ c n ] of modes sat-urate at late times, leading to permanent propagatingprofiles with ξ nk ( t ) = c n e inω ( k ) t (see below).Next, we perform fully nonlinear simulations to vali-date the weakly nonlinear predictions. The vortex sheetmethod is a standard sharp-interface technique for sim-ulating dynamics of Hele-Shaw flows. It is based on aboundary integral formulation in which the interface isformally replaced by a generalized vortex sheet [21] witha distribution of vortex strengths γ ( s, t ), where s is thearclength coordinate. We adapt this approach to handleferrofluids under imposed magnetic fields. First, we ex-press the velocity of the interface solely in terms of theinterface position. To do so, it is convenient to identifythe position vector in R with a scalar z ( s, t ) ∈ C ( ∗ de-notes complex conjugate) [21–23]. Second, to advancethe interface, we solve the dimensionless equations z ∗ t = − γ z s + 12 πi P (cid:73) γ ( s (cid:48) , t ) z ( s, t ) − z ( s (cid:48) , t ) ds (cid:48) , (7a) γ (cid:60) (cid:26) z s πi P (cid:73) γ ( s (cid:48) , t ) z ( s, t ) − z ( s (cid:48) , t ) ds (cid:48) (cid:27) + (cid:2) κ ( s, t ) − ( M · ˆn ) − Ψ (cid:3) s , (7b)iteratively for the velocity z t , where ( · ) t ≡ ∂ ( · ) /∂t ,( · ) s ≡ ∂ ( · ) /∂s , i = √−
1, and P represents principalvalue integration. Here, Ψ = N Ba r − + N Br r , and( M · ˆn ) = χ (cid:2) √ N Ba r − ( ˆe θ · ˆn ) + √ N Br r ( ˆe r · ˆn ) (cid:3) is thedimensionless magnetic normal stress. Time advance-ment is accomplished by the Crank–Nicolson scheme.The spatial discretization is implemented on an array ofLagrangian points ( N = 1024) with uniform ∆ s [24]. Evolutionary dynamics
The evolution of perturbedharmonic modes ξ k under the fully nonlinear simulationand the weakly nonlinear approximation are shown inFig. 2(a). Starting from small initial values ( ξ k | k =7 =0 . ξ nk = 0 for n >
1) with N Ba , N Br , R, χ set sothat the most unstable mode is equal to the fundamentalmode ( k m = k = 7), they saturate at late times. The per-turbed circular interface grows exponentially in the linearregime and then matches the weakly nonlinear approxi-mation at intermediate times ( t ∈ [0 , t w ]). The nonlinearsimulations take longer to saturate ( t ∈ [ t w , t e ]) and doso at higher final amplitudes compared to the weak non-linear result. The time-domain evolution is also shown inFig. 2(b), evolving from a nearly flat (unwound circular)interface into a permanent propagating profile [25].The rotating droplet, shown in Fig. 2(c), has a polyg-onal shape with the symmetry set by the fundamentalmode, k = 7. The fully nonlinear profile has a sharperpeak compared to the weakly nonlinear approximation,which is otherwise in good agreement. The key discoveryof the present work is the stable rotating shape , which wenow seek to analyze as a nonlinear wave phenomenon. When does weakly nonlinear stability imply nonlinearstability?
One of deficiencies of linear and weakly non-linear analyses is that they do not provide sufficient con-ditions for stability. Linearly stable base states can benonlinearly unstable [26], and vice versa . Importantly,however, our nonlinear traveling wave solution is a local attractor (following the terminology from [27]), as shownin Fig. 2(d).Shapes in a neighborhood of the propagating pro-file, subject to small ( ξ k, k /ξ fk, k (cid:28)
1) or intermediate( ξ k, k /ξ fk, k = O (1)) initial perturbations, converge toit. Larger perturbations (shaded region) lead to non-linear instability of the weakly nonlinearly stable pro-files; infinitely long rotating “fingers” form. Convergenceto the attractor is sensitive to the initial amplitude ofthe first harmonic mode ξ k . For the chosen parameters, λ ( k ) > λ (2 k ) <
0: high wavenumbers decay andthe fundamental wavenumber grow in the linear stage.Consequently, for low ξ k /ξ fk and high ξ k /ξ f k , the lowwavenumber modes grow and saturate, as high wavenum-ber modes decay exponentially in the linear regime. Withhigher initial ξ k /ξ fk , the perturbed droplet will not go FIG. 2. (a) The evolution of the first 5 harmonic modes fromfully nonlinear simulation (solid) and weakly nonlinear ap-proximation (dashed), for N Ba = 1 .
0, N Br = 37, R = 1, and χ = 1 (same parameters for (b), (c) and (d)). (b) The fullynonlinear evolution of the interface from a small perturba-tion of the flat base state into a permanent traveling wave(rotating droplet). (c) Comparison between the final shapefrom fully nonlinear simulation (solid) and weakly nonlinearapproximation (dashed); (d) Stability diagram based on thefirst two harmonic modes of the final shape (marked with (cid:78) )shown in (b); ◦ (resp. × ) denotes the stable (resp. unstable)initial conditions, the unstable region is shaded, and the ‘ f ’superscript represents the final harmonic mode amplitude. through the linear region, and the amplitudes of bothmodes will rapidly increase to create a skewed shape,with multivalued h ( θ, t ) for which harmonic modes canno longer be defined. Note that Fig. 2(d) is a projec-tion in the ( ξ k , ξ k ) plane, where the initial values of ξ k , ξ k , ξ k are set as the final amplitudes (and phases)from the weakly nonlinear equations [28]. Propagation velocity
A permanent traveling waveprofile has ξ ( θ, t ) = Ξ( kθ − ωt ), and v f = ω/k isits propagation velocity. Expressing the modes’ com-plex amplitudes as ξ nk ( t ) = c n e − inω ( k ) t , with constant c n ∈ C that account for their relative phases, we have v p ( k, t ) = nω ( k ) /nk = v f . The mean v p of the first fiveharmonics is used to calculate v Ff for the fully nonlinearsimulation and also v Wf for the weakly nonlinear approx-imation. Meanwhile, v Lf = v p as given below Eq. (6).For a quantitative comparison, three sets of parametersare considered, fixing χ = 1. Two sets (i) and (ii) are for k m = 7, and the variation of N Br is according to Eq. (6).A third set (iii) explores the effect of k m under the samelinear propagation velocity v Lf . Figure 3(a) compares thefinal propagating velocity predictions. Both v Lf and v Wf are in relatively good agreement with v Ff for small ve-locities. When N Ba → v f = 0)exist [12]. For higher v f , the larger deviation in the pre-dictions highlights the importance of nonlinearity. Nev-ertheless, the linear and weakly nonlinear results follow asimilar trend. Importantly, v Lf and v Wf help identify thekey control factors: the coupled magnetic field strength √ N Ba N Br and the radius of the initial droplet R . Thesalient physics uncovered is that the propagating veloc-ity can be non-invasively tuned. Traveling wave shape
The most unstable mode k m sets the propagating profile, which has a sharper peakfor higher k m [Fig. 3(c)]. To quantify the shape change,we introduce the skewness Sk ( t ) = (cid:104) ξ (cid:105) / (cid:104) ξ (cid:105) / , which isused to define the vertical asymmetry of nonlinear sur-face water waves [29, 30]; Sk > (cid:104) · (cid:105) = π (cid:82) π ( · ) dθ . Fig-ure 3(b) shows that Sk (for the fully nonlinear propa-gation) increases with k m , as expected from the sharperpeaks in Fig. 3(c). This observation also explains why v Lf becomes a worse approximation of v Ff as k m increases[inset of Fig. 3(a)]: smoother peaks (lower k m ) are bettercaptured by the linear theory based on harmonic modes.Figure 3(c) reveals that the wave profile for k m = 5(N Ba = 1 .
9, N Br = 19 .
5) is more asymmetric than theone for k m = 9 (N Ba = 0 .
60, N Br = 60 . Ba = 0), the stationary shape hasazimuthal symmetry [12]. To understand the asymmetryof propagating shapes induced by the combined magneticfield [31], we extend the parameters of case (i) to a newcase (iv): k m = 7, R = 1 and N Br varying according toN Ba (see Fig. 4 caption). To quantify the fore-aft asym-metry of the shape, we introduce As ( t ) = (cid:104)H [ ξ ] (cid:105) / (cid:104) ξ (cid:105) / [29, 30]; H [ · ] is the Hilbert transform. For As >
0, wavestilt “forward” ( i.e. , counter-clockwise).Figure 4(a) shows As ( t ) for different √ N Ba N Br , whichquantifies the coupled field effect, starting with smallsymmetric perturbations. For a stable case, As ( t ) reachesa maximum value ( t ≈ t ) during the initial unstableweak nonlinear growth (dark shadow region), and asymp-totes to a value close to zero ( t ≥ t ). The differences inthe final propagating profile (under the same k m ) shownin Fig. 4(b) are hard to capture, which is consistent withthe observation in Fig. 3(b). For the unstable cases,“wave breaking” occurs, which is highlighted by a changeof sign of As . Also, now, As ( t ) no longer saturates at late t . Instead As ( t ) crosses zero (at t (cid:38) t ) and approaches a FIG. 3. (a) Comparison of the propagation velocity predictedby linear theory v Lf (dashed), weakly nonlinear theory v Wf (empty symbols), and fully nonlinear simulation v Ff (filledsymbols). The circles represent results for case (i) R = 1 fixedand N Ba ∈ [0 , . , . , . , , , Ba = 1 fixed with N Br varying accordingto R ∈ [0 . , . , . , . k m ∈ [5 , , , , R = 1 and N Ba , N Br determined so that v Lf = 85 .
16. (b) The skewness Sk of the fully nonlinear profile.(c) The permanent wave shape (only one wavelength shown). singularity. This unstable example is shown in Fig. 4(c)and the second column of Fig. 4(d). As its amplitudefirst grows, the wave tilts forward ( t = t ), but nonlineareffects restore its symmetry ( t = t ). Subsequently, thewave tilts backwards ( t = t , t ) and its amplitude con-tinues to grow ( t = t ) [32]. The distorted wave has awider base and evolves into long unstable fingers.Note that N Ba / N Br also increases with √ N Ba N Br forour choices of N Ba and N Br . Equation (4) shows thatthe radial magnetic field is destabilizing, while surfacetension ( k > k m , increasing N Ba / N Br can induce instabilitybecause it engenders a larger v f (and As ), leading to aglobal bifurcation with Fig. 2(d) as one stable slice. Thisresult has an analogy to solitary waves in equations of theKortweg–de Vries (KdV) type. Specifically, initial per-turbations grow, deforming a shape until nonlinearity isbalanced by dispersion, when a permanent wave emerges[33]. However, depending on the form of the nonlinearity,not all such permanent waves are stable attractors, andconditions must be placed on the wave speed [34]. Discussion
Although the manipulation of the lineargrowth rate of interfacial perturbations in Hele-Shawcells is well studied [35], the control of the dynamic nonlinear patterns is not. Our approach harnesses themagnitude and the direction of coupled magnetic fieldsto generate ferrofluid droplets, with well-characterized
FIG. 4. (a) Time evolution of the wave profile asymme-try for different combination of N Ba ∈ [0 , . , , , , , , , Br varying so that k m = 7 for R = 1. Solid curvesrepresent stable cases yielding a propagating profile; dashedcurves represent unstable cases in which the profile distortsand grows without bound. (b) Asymptotic stable propagat-ing wave profile. (c) The growing unstable interfacial pat-tern (N Ba = 8 , N Br = 41 at t = t ) (d) Stable (left, withN Ba = 1 , N Br = 37) and unstable (right, N Ba = 8 , N Br = 41)evolution of the profile. Corresponding skewness evolutioncurve and the time are marked with white dot on (a). shapes and rotational speeds, by purely external means.Open questions remain: e.g. , which fundamentalmodes evolve into propagating shapes? Work on the sta-tionary problem [15, 36] gives a hint, however, for a prop-agating shape the Birkhoff integral equation [37] must besolved, making an extension of [15, 36] challenging. In-terestingly, our simulations also reveal that patterns pre-dicted as stable by weakly nonlinear analysis can be un-stable [38]. Additionally, does this system accommodatemore than one propagating wave? If so, do such waveskeep their shapes upon collision, as with soliton interac-tions [33, 39]? Previous studies derived KdV equationsfor unidirectional small-amplitude, long-wavelength dis-turbances on fluid–fluid interfaces in Hele-Shaw [40] andaxisymmetric ferrofluid configurations [41, 42], demon-strating the celebrated “sech ” solitary wave. Instead,in our study without such restrictions, we discovered pe-riodic traveling nonlinear waves, which are akin to the cnoidal solutions of periodic KdV, i.e. , the fundamentalnonlinear modes (“soliton basis states”) [43]. Addition-ally, we observed wave breaking [Fig. 4(d)].Finally, it would be of interest to verify the pro-posed shape manipulation strategies by laboratory ex-periments. Previous theoretical studies [15, 44–46] sug-gest that many exact stationary droplet shapes are un-stable, thus their relevance to experimental studies is lim-ited. On the other hand, the nonlinear simulations in ourstudy, showing stable rotation, pave the way for futureexperimental realizations. Acknowledgements
This material is based upon worksupported by the National Science Foundation underGrant No. CMMI-2029540. ∗ [email protected] † [email protected][1] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B.Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrody-namics of soft active matter, Rev. Mod. Phys. , 1143(2013).[2] D. Saintillan, Rheology of Active Fluids, Annu. Rev.Fluid Mech. , 563 (2018).[3] A. Sokolov, M. M. Apodaca, B. A. Grzybowski, and I. S.Aranson, Swimming bacteria power microscopic gears,Proc. Natl Acad. Sci. USA , 969 (2010).[4] G. Z. Lum, Z. Ye, X. Dong, H. Marvi, O. Erin, W. Hu,and M. Sitti, Shape-programmable magnetic soft matter,Proc. Natl Acad. Sci. USA , E6007 (2016).[5] T. Xu, J. Zhang, M. Salehizadeh, O. Onaizah, andE. Diller, Millimeter-scale flexible robots with pro-grammable three-dimensional magnetization and mo-tions, Science Robotics , eaav4494 (2019).[6] J. Giltinan, P. Katsamba, W. Wang, E. Lauga, andM. Sitti, Selectively controlled magnetic microrobotswith opposing helices, Appl. Phys. Lett. , 134101(2020).[7] B. J. Nelson, I. K. Kaliakatsos, and J. J. Abbott, Mi-crorobots for Minimally Invasive Medicine, Annu. Rev.Biomed. Eng. , 55 (2010).[8] R. Rosensweig, Ferrohydrodynamics (Dover Publications,Mineola, NY, 2014) republication of the 1997 edition.[9] E. Blums, A. Cebers, and M. Maiorov,
Magnetic Fluids (De Gruyter, Berlin, 1997).[10] R. E. Rosensweig, M. Zahn, and R. Shumovich,Labyrinthine instability in magnetic and dielectric fluids,J. Magn. Magn. Mater. , 127 (1983).[11] D. P. Jackson, R. E. Goldstein, and A. O. Cebers, Hydro-dynamics of fingering instabilities in dipolar fluids, Phys.Rev. E , 298 (1994).[12] P. H. A. Anjos, S. A. Lira, and J. A. Miranda, Fingeringpatterns in magnetic fluids: Perturbative solutions andthe stability of exact stationary shapes, Phys. Rev. Fluids , 044002 (2018).[13] A. L. Bertozzi and M. C. Pugh, Long-wave instabilitiesand saturation in thin film equations, Comm. Pure Appl.Math. , 625 (1998).[14] D. Bensimon, L. P. Kadanoff, S. Liang, B. I. Shraiman,and C. Tang, Viscous flows in two dimensions, Rev. Mod.Phys. , 977 (1986).[15] S. A. Lira and J. A. Miranda, Ferrofluid patterns in Hele-Shaw cells: Exact, stable, stationary shape solutions,Phys. Rev. E , 013129 (2016).[16] J. A. Miranda and R. M. Oliveira, Time-dependent gapHele-Shaw cell with a ferrofluid: Evidence for an interfa-cial singularity inhibition by a magnetic field, Phys. Rev.E , 066312 (2004).[17] So, ∇ | H | (cid:54) = is the main contribution to the body force,and the demagnetizing field is negligible [15, 16, 41, 47].[18] J. A. Miranda and M. Widom, Radial fingering in a Hele-Shaw cell: a weakly nonlinear analysis, Physica D ,315 (1998).[19] See Supplemental Material at [URL will be inserted bypublisher] for the expressions of the functions F and G .[20] Intuitively, from vector projection, we observe that onlythe combined azimuthal and radial magnetic field can break the symmetry and cause a force imbalance leadingto motion.[21] A. Prosperetti, Boundary Integral Methods, in Drop-Surface Interactions , CISM International Centre for Me-chanical Sciences, Vol. 456, edited by M. Rein (Springer-Verlag, Wien, 2002) pp. 219–235.[22] M. J. Shelley, A study of singularity formation in vortex-sheet motion by a spectrally accurate vortex method, J.Fluid Mech. , 493 (1992).[23] M. J. Shelley, F.-R. Tian, and K. Wlodarski, Hele-Shawflow and pattern formation in a time-dependent gap,Nonlinearity , 1471 (1997).[24] See Supplemental Material at [URL will be insertedby publisher] for further details, including algorithmflowchart and grid convergence study.[25] See Supplemental Material at [URL will be inserted bypublisher] for a video of this process.[26] P. G. Kevrekidis, D. E. Pelinovsky, and A. Saxena, WhenLinear Stability Does Not Exclude Nonlinear Instability,Phys. Rev. Lett. , 214101 (2015).[27] S. Li, J. S. Lowengrub, J. Fontana, and P. Palffy-Muhoray, Control of Viscous Fingering Patterns in a Ra-dial Hele-Shaw Cell, Phys. Rev. Lett. , 174501 (2009).[28] Even though Fig. 2(a) indicates ξ k makes a non-trivialcontribution to the final shape, while ξ k , ξ k play asmaller role, the projection is sufficient to conclude thatthe propagating wave profile is an attractor.[29] A. B. Kennedy, Q. Chen, J. T. Kirby, and R. A. Dalrym-ple, Boussinesq modeling of wave transformation, break-ing, and runup. I: 1D, J. Waterw. Port. Coast. , 39(2000).[30] T. J. Maccarone, The biphase explained: understandingthe asymmetries in coupled Fourier components of astro-nomical time series, Mon. Not. R. Astron. Soc. , 3547(2013).[31] See Supplemental Material at [URL will be inserted bypublisher] for a comment on the source of the shapeasymmetry on the basis of the magnetic normal stressexpression.[32] The calculation of As now fails because H requires theperturbation ξ ( θ, t ) to be single-valued in θ .[33] N. J. Zabusky and M. D. Kruskal, Interaction of “Soli-tons” in a Collisionless Plasma and the Recurrence ofInitial States, Phys. Rev. Lett. , 240 (1965).[34] J. L. Bona, P. E. Souganidis, and W. A. Strauss, Stabil-ity and instability of solitary waves of Korteweg-de Vries type, Proc. R. Soc. Lond. A , 395 (1987).[35] L. C. Morrow, T. J. Moroney, and S. W. McCue, Numer-ical investigation of controlling interfacial instabilities innon-standard Hele-Shaw configurations, J. Fluid Mech. , 1063 (2019).[36] E. S. G. Leandro, R. M. Oliveira, and J. A. Miranda, Ge-ometric approach to stationary shapes in rotating Hele-Shaw flows, Physica D , 652 (2008).[37] G. Birkhoff, Taylor instability and laminar mixing , Tech.Rep. LA-1862 (Los Alamos Scientific Laboratory, 1954).[38] See Supplemental Material at [URL will be inserted bypublisher] for an example showing that perturbationswith k = 4 will not evolve into either a stationary or apropagating shape (although both are predicted to existby weakly nonlinear analysis).[39] T. Soomere, Solitons Interactions, in Encyclopedia ofComplexity and Systems Science , edited by R. A. Meyers(Springer, New York, NY, 2009) pp. 8479–8504.[40] M. Zeybek and Y. C. Yortsos, Long waves in parallel flowin Hele-Shaw cells, Phys. Rev. Lett. , 1430 (1991).[41] D. Rannacher and A. Engel, Cylindrical KortewegdeVries solitons on a ferrofluid surface, New. J. Phys. ,108 (2006).[42] E. Bourdin, J.-C. Bacri, and E. Falcon, Observation ofAxisymmetric Solitary Waves on the Surface of a Fer-rofluid, Phys. Rev. Lett. , 094502 (2010).[43] A. R. Osborne, E. Segre, G. Boffetta, and L. Cavaleri,Soliton basis states in shallow-water ocean surface waves,Phys. Rev. Lett. , 592 (1991).[44] S. A. Lira, J. A. Miranda, and R. M. Oliveira, Station-ary shapes of confined rotating magnetic liquid droplets,Phys. Rev. E , 036318 (2010).[45] E. O. Dias, S. A. Lira, and J. A. Miranda, Interfacialpatterns in magnetorheological fluids: Azimuthal field-induced structures, Phys. Rev. E , 023003 (2015).[46] E. ´Alvarez-Lacalle, J. Ort´ın, and J. Casademunt, Non-linear Saffman-Taylor instability, Phys. Rev. Lett. ,054501 (2004).[47] P. H. A. Anjos, G. D. Carvalho, S. A. Lira, and J. A.Miranda, Wrinkling and folding patterns in a confinedferrofluid droplet with an elastic interface, Phys. Rev. E , 022608 (2019). Supplementary Material forTuning a Magnetic Field to Generate Spinning Ferrofluid Droplets with ControllableSpeed via Nonlinear Periodic Interfacial Waves
Zongxin Yu, ∗ and Ivan C. Christov, † School of Mechanical Engineering, Purdue University, West Lafayette, Indiana 47907, USA ∗ Electronic address: [email protected] † Electronic address: [email protected](Dated: September 11, 2020)
SECOND-ORDER MODE-COUPLING TERMS
The functions describing the second-order mode-coupling terms in Eq. (3) in the main text are found to be: F ( k, k (cid:48) ) = | k | R (cid:26) N Ba R [3 − χk (cid:48) ( k − k (cid:48) )] + N Br { χ [ k (cid:48) ( k − k (cid:48) ) + 1] }− R (cid:20) − k (cid:48) k (cid:48) + k ) (cid:21) + 2 χ √ N Ba N Br R ik (cid:48) (cid:27) , (S1a) G ( k, k (cid:48) ) = 1 R [(sgn( kk (cid:48) ) − | k | − , (S1b)where sgn( x ) = x/ | x | for x (cid:54) = 0 and sgn(0) = 0. SOURCE OF THE SHAPE ASYMMETRY
The dimensionless governing Eq. (1) and pressure boundary condition in Eq. (2) can be rewritten as v = − ∇ (cid:18) p − N Ba r − N Br r (cid:19) , and (S2) p = κ − (cid:20) χ N Ba r ( ˆe θ · ˆn ) + χ N Br r ( ˆe r · ˆn ) + 2 χ (cid:112) N Ba N Br ( ˆe θ · ˆn )( ˆe r · ˆn ) (cid:21) , (S3)where ˆe θ · ˆn = − h θ ( h + h θ ), ˆe r · ˆn = h/ ( h + h θ ). The magnetic scalar potential in Eq. (S2) results from thebody force, and the terms pre-multiplied by χ in Eq. (S3) represent the magnetic normal stress. For a droplet withsymmetric azimuthal perturbation, the body force alone cannot break the symmetry. Therefore, the asymmetry ofshapes discussed in the main text is to be attributed to the magnetic normal stress.This observation can be intuitively understood by considering one wavelength of a symmetric waveform. The firstthree terms on the right-hand side of Eq. (S3) are equal on both sides of the peak, while the fourth term changesat the peak due to the sign of h θ , which requires different curvatures on either side of the peak to remain balanced.Therefore, √ N Ba N Br can be taken as the measure of the coupling effect between the magnetic field components. UNSTABLE PATTERN WITH k = 4 As mentioned in the main text, linear (or even weakly nonlinear) stability does not always imply nonlinear stability(see also [26]). Indeed, it is not known under what conditions the weakly-nonlinear stable droplet shapes are actuallynonlinearly stable. In the main text, we presented examples for which this implication holds true. Here, in Fig. S1,we demonstrate, for completeness, an example to the contrary. To the best of our knowledge, such an example hasnot been analyzed before, and thus remains an avenue of future work. This exploration also requires caution, to ruleout physical from numerical instability.
FIG. S1. The unstable evolution of the first 5 harmonic modes ( k = 4 , , , ,
20) from fully nonlinear simulation (solid) andtheir stable evolution from weakly nonlinear approximation (dashed) for (a) nonrotating ( k m = 4, N Ba = 0) and (b) rotating( k m = 4, N Ba = 1) shapes. IMPLEMENTATION OF THE NUMERICAL METHOD AND GRID CONVERGENCE
The principal value integration in Eqs. (7) is performed numerically by a spectrally accurate spatial scheme [22]:
P V j = 2∆ s πi (cid:88) j + k odd γ k z j − z k , (S4)where a j subscript denotes the evaluation of a quantity at the j th Lagrangian grid point j ∆ s with ∆ s = L/N , L = (cid:72) ds , and N is the number of grid points. The parametrization of the interface via its arclength reduces thestiffness of the numerical problem caused by the presence of third-order spatial derivatives. A rearrangement of thegrid points is conducted with cubic interpolation, after each time step, to maintain uniform grid spacing ∆ s . Theuniform arclength spacing then allows the use of the second-order central differentiation formulæ for all derivatives.A fixed-point iteration scheme is used to resolve the implicit Eq. (7b) to obtain γ j at each interface point z j , as shownschematically in Fig. S2.Time advancement (superscripts denote the time step number) is accomplished with a Crank–Nicolson scheme: z ∗ n +1 = z ∗ n − ∆ t (cid:20) γ n +1 z n +1 s + γ n z ns (cid:21) + ∆ t (cid:2) P V n +1 + P V n (cid:3) , (S5)where both of the nonlinear terms γ n +1 / z n +1 s and P V n +1 are obtained by subiteration with index m as shown inFig. S2. Equation (S5) converges and z n +1 = z m +1 when (cid:107) z m +1 − z m (cid:107) < tol z with tol z = 0 . j | z nt | ∆ t .A grid convergence study with 4 levels of the grid resolution was conducted for two cases: k m = 7 and k m = 9 withN Ba = 1. The most frequently used case in the main text is k m = 7, while the sharper peaks for k m = 9 demandon the highest grid resolution. Figure S3(a) shows the spectral energy content of harmonic modes ( k , 2 k , 3 k , . . . )of the propagating waveform, where the “piling up” near the tail on the finest grid ( N = 2048) is numerical noise.This plot supports our decision to consider the N = 1024 grid as offering sufficient resolution. Figure S3(b) showsthe root-mean-squared error in the shape z itself, taking ˆ z as the “reference shape” on the N = 2048 grid. The errordecreases with grid refinement. The error at N = 256 for k m = 9 is not shown for the propagating shape because thescheme is not even stable on such a coarse mesh for this case.Figure S3 further shows the evolution of (c) the skewness Sk and (d) the asymmetry As . The skewness matcheswell on all grids used, showing it is a well-converged quantity, while the asymmetry is seen to be more sensitive to thegrid resolution. The differences between N = 1024 and N = 2048 are small enough so that it is safe to use N = 1024for the simulations reported in the main text, considering the significantly higher computational cost incurred byusing finer grids. FIG. S2. Flow chart of the vortex-sheet algorithm using Crank–Nicolson for time advancement and fixed-point iteration forresolving the implicit nonlinear terms.FIG. S3. Grid convergence study for fundamental modes k m = 7 (black) and k m = 9 (blue) with N Ba = 1 and N = 256(dotted), N = 512 (dot-dashed), N = 1024 (dashed), and N = 2048 (solid). (a) Spectral energy of harmonic modes ( k , 2 k , 3 k , . . . ). (b) The root-mean-square error taking the N = 2048 solution ˆ z as “exact”. (c) Grid convergence of the evolution of theskewness Sk ( t ). (d) Grid convergence of the evolution of the asymmetry As ( tt