Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
Open AccessResearch articles

Stability of an inhomogeneous ferrofluid in a channel, subject to a normal field

S. H. Ferguson Briggs

S. H. Ferguson Briggs

Department of Mathematics, Imperial College London, Exhibition Road, South Kensington, London SW7 2BX, UK

[email protected]

Contribution: Formal analysis, Investigation, Resources, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

A. J. Mestel

A. J. Mestel

Department of Mathematics, Imperial College London, Exhibition Road, South Kensington, London SW7 2BX, UK

Contribution: Formal analysis, Investigation, Project administration, Supervision, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed



    The stability of a ferrofluid with a fairly arbitrary non-uniform magnetic susceptibility between two parallel walls, subject to a non-uniform magnetic field acting normal to the walls, is investigated. The susceptibility may depend on position and the field strength, and the stationary state is such that the gradient of the susceptibility with respect to the modulus of the field is negative. Previous work suggests that the configuration may be linearly unstable, as regions of higher susceptibility do not coincide with regions of strongest field, and this is proved here. Adding a constant field in the plane of the layer suppresses parallel instabilities, but has no effect on those orthogonal to it. However, it is demonstrated that stability can be achieved by applying a rapidly rotating field.

    1. Introduction

    Ferrofluids (FFs) are colloidal, nano-size magnetic particles suspended in a carrier liquid [1], where the particles are stable against sedimentation and agglomeration [2]. The initial magnetic susceptibility of an FF is 104 times larger than an ordinary paramagnetic fluid [2,3], resulting in them being extremely receptive to an imposed magnetic field, which is useful in many disciplines. Most common are their industrial applications for dynamic sealing and heat dissipation, for example, in audio speakers [4]. Bio-medical applications include magnetic drug targeting, where the FF allows for a localized treatment, and hyperthermia treatment in cancer therapy [58]. Theoretical analysis of FFs is, therefore, useful in determining the stability of various geometrical configurations and the possible parameter ranges for a chosen application. Previous models of FF systems focus on FFs with homogeneous magnetic susceptibilities, but in the present paper we consider FFs with inhomogeneous susceptibilities, for which equilibria do exist, in an effort to produce a more general model relevant for current applications and to influence new applications.

    We investigate the linear stability of an FF between two planar walls, subject to a magnetic field applied normal to the walls. Following the work by Rosensweig [1], we assume a co-linear relationship between the magnetization M of the FF and the applied field H, such that M=χH, now named the quasi-equilibrium approach [9], where χ is the magnetic susceptibility. We allow χ to depend on both the magnitude of the field H and position, and could, therefore, be applied to a nonlinear susceptibility and spatial variation. We are particularly interested in configurations when the magnetic field strength and susceptibility are functions of distance from the walls. The simplest example consists of an FF layer with a non-magnetizable fluid above it. It is well known that the interface between the two fluids is destabilized by a normal magnetic field [10]. In an analogous electro-hydrodynamic (EHD) problem, Li et al. [11] investigate the stability of an interface between two dielectric fluid layers between two walls. They find the configuration is unstable, a result that holds for an analogous system of two FFs with constant magnetic susceptibilities. The direct analogue between dielectrics and FFs is well known for ferrofluids exhibiting a linear susceptibility [1,12]. Although the analogy holds for nonlinear susceptibilities too [13], the EHD literature for a dielectric exhibiting nonlinear polarization is limited. Yecko [14,15] consider an analogous set up to Li et al. [11] but with the lower fluid a FF with nonlinear susceptibility, and apply the Langevin function for paramagnetic systems to describe the magnetization of the FF. For a nonlinear (and a non-constant) susceptibility, the magnetic forcing is felt in the bulk of the fluids, but Yecko [14,15] assume a form for the magnetic stress tensor such that the forcing remains confined to the interface to show the system is unstable. Zelazo & Melcher [12] consider the interfacial instability of two FFs with nonlinear susceptibility between two current sheets, subject to a spatially varying field in a two-dimensional cylindrical domain. They implement a quasi-one-dimensional model, transforming a cylindrical geometry to a planar geometry, and use a Boussinesq approximation to obtain a dispersion relation, showing that when the field decreases upwards, and the stronger FF is the lower fluid, only capillary forces are destabilizing. Yet, if the stronger FF is the upper fluid, both capillary and magnetic forcing are destabilizing, implying that the magnetic force is destabilizing when the highest region of magnetic susceptibility does not coincide with the highest region offield.

    Ferguson Briggs & Mestel [16] (henceforth FBM) argue that given a stationary state with a non-uniform field and susceptibility, such that the gradient of the susceptibility of the fluid with respect to the gradient of the magnitude of the field, dχ/dH, is negative, the system may be unstable. This was proven for a cylindrical domain where the susceptibility of a column of FF varied continuously with radius, subject to an azimuthal field which decreased as the reciprocal of the radius, driven by an axial current. The configuration was shown to be unstable if χ increased radially at any point in the fluid, but the instability could be suppressed by the addition of another field. The need for a central wire carrying the axial current is inconvenient for many applications. This motivates the study of the stability of equilibria with non-uniform susceptibility and field strength in other geometries.

    In this paper, we consider a channel with a fairly arbitrary one-dimensional distribution of magnetic particles giving rise to the susceptibility χ(z), where z is normal to the channel walls. If this FF is placed in a uniform magnetic field in the z-direction, the field will adjust to an equilibrium strength H(z), for which it transpires that dχ/dH is negative. Hence from FBM and the conclusions of Zelazo & Melcher [12], we postulate that the system would be unstable. In §3, we prove the system is indeed linearly unstable to all modes and in §4 we seek to stabilize the system. It is known that a field applied parallel to the plane interface between an FF and a non-magnetic fluid will suppress linearly unstable modes at the interface whose wavevector is parallel to the field. Yet, for a three-dimensional system, modes perpendicular to the field will not be dampened [17]. Korovin [18] considers the stability in the inviscid limit of the classical Rosensweig instability, subject to a tilted piecewise-constant magnetic field, and shows that in the presence of a horizontal field, a stronger vertical field is needed to produce the instability. Dorbolo & Falcon [19] show experimentally (and theoretically) a horizontal magnetic field acting on sinusoidal waves at a fixed frequency produces a monotonous dispersion relation, but observe that in the nonlinear regime wave turbulence occurs. In a channel system, Zelazo & Melcher [12] and Yecko [15] both show adding a constant field down the channel does not stabilize all modes, and we demonstrate that a constant field across the channel is not sufficient in stabilizing the system outlined here either. Schumacher et al. [20] analyse the nonlinear effects from applying a steady or oscillating field in the direction of a pressure-driven flow in a three-dimensional, wall-bounded FF system. They find the oscillating field requires less of a pressure drop increase to maintain a constant flow rate than a steady field. Rannacher & Engel [17] employ a rotating field to stabilize the classical Rayleigh–Taylor instability, for an FF as the upper (more dense) fluid. By Floquet theory, they find a rotating field will stabilize the modes where the modulus of the wavenumber satisfies a threshold value, but otherwise the modes remain unstable. However, we prove a sufficiently strong, rapidly rotating field does stabilize all modes for our system in §4. We, therefore, have a mechanism for maintaining an arbitrary layered distribution of magnetic particles using a non-invasive external magnetic field.

    2. Magnetic force and stress tensor

    We assume the ferrofluid is an electric insulator and that we have an isothermal system throughout. We neglect particle interactions and magnetoviscous effects, whereby the rotation of individual particles when the vorticity and field are not aligned could lead to variation in the fluid viscosity [3,21]. We employ a quasi-equilibrium approach, where the magnetization M is attained instantaneously and is parallel to the field H. Consequently, M=χH, where χ is not necessarily constant. The field satisfies

    and the induced field, B, satisfies
    Equation (2.1) allows us to define a magnetic potential, ϕ, such that
    Due to co-linearity, B=μ0(1+χ)H, where μ0 is the magnetic permeability of free space, and therefore
    We require continuity of the normal component of B and the tangential component of H through the walls.

    The stress tensor for a Newtonian FF, neglecting magneto-strictive effects on the flow, is given by

    where H=|H|, η is the viscosity, p is the pressure and u is the velocity of the fluid.The force density for a FF subject to H is
    and if χ does not depend explicitly on H,
    f appears in the Navier–Stokes equation as
    where ρ is the density of the fluid [1]. We assume the fluid is incompressible such that
    and taking the curl of (2.8) gives the vorticity equation
    where ω=×u is the vorticity of the flow. Equation (2.10) holds for both (2.6) and (2.7) and it follows that for a stationary state we require χχ(H). Note that χχ(H) can hold when χ and H both depend on position, and they do not need to depend on one another explicitly. Similarly to Zelazo & Melcher [12], and as in FBM we assume,
    such that a displaced fluid parcel retains its value of χ and its physical properties over time scales of interest. If χ depends explicitly on H, then (2.11) allows for an initial dependence of χ on H, and assumes that any subsequent change in χ due to its dependence on H, will occur on a slower time scale than times scales of interest.

    3. Formulation

    Stationary equilibria must satisfy both χχ(H) and (2.4). The simplest solutions have χ constant throughout the fluid. Other solutions are hard to find, but a more complex solution family of (2.10) has

    for a positive constant A, such that A/H1 throughout the fluid, ensuring χ0. χ does not need to depend explicitly on H for (3.1) to be satisfied. It should be noted that (3.1) is not a physical or chemical relationship between H and χ, and is not valid for all H (for example H0 and H). Rather, it is a formula which happens to hold for a particular equilibrium. If χ takes an arbitrary but specific distribution in space, a field applied to the FF will adapt spatially to satisfy (2.2) and (2.8). For χ given by (3.1), (2.4) becomes
    Since ϕ/|ϕ| is a unit normal to the surface of constant ϕ, (3.2) implies that the equipotential surfaces have zero mean curvature and hence minimize their area [22]. Consequently, the governing equations (subject to boundary conditions) are satisfied for a stationary fluid when the spatial distributions of χ and H satisfy (3.1), and equipotential surfaces have zero mean curvature.

    The simplest zero-curvature surfaces are planar and in this paper we investigate the one-dimensional equilibria with H=H0(z)ez, where χ is given by (3.1). We consider a three-dimensional, cartesian coordinate system, where an FF of spatially varying χ is between two solid plates at z=z1 and z=z2, where z points normal to the plates and x,y perpendicular to each other along the plates. The FF is incompressible and isothermal with constant ρ, η and μ0. The fluid is initially at rest, and the susceptibility of the fluid satisfies

    but elsewhere,
    where A/H01. It follows that
    where l=1,2 for zz1,zz2, respectively. A schematic of the problem is shown in figure 1. Here, we consider rigid boundaries, but we could have no plates and impose the appropriate boundary conditions as z±.
    Figure 1.

    Figure 1. Schematic of the system. (Online version in colour.)

    Interestingly, given (3.1), it follows that


    As discussed in FBM, such equilibria may release magnetic energy by re-distributing their ferrous particles and so may be prone to instability. We illustrate this argument in figure 2. Figure 2a shows the parcels of fluid with highest susceptibility do not coincide with the regions of strongest field. Two parcels are swapped, resulting in the location of each parcel, and its corresponding value of susceptibility, ‘coinciding more’ with the strength of the field, shown in figure 2b. This may result in a release of energy, as less energy is required for this formation, and the perturbation can use this energy to grow. The lowest energy formation is shown in figure 2c.

    Figure 2.

    Figure 2. An illustration of an energy argument: the field acts in the vertical direction, decreasing in strength with increasing height, where the darker the blue colour, the stronger the field. The FF is illustrated as parcels with individual values of susceptibility, where the darker the orange colour, the higher the susceptibility. (a) shows the upper parcel has the highest susceptibility, the bottom parcel has the least, and the intermediate parcels are such that the susceptibility increases with height. Energy is released moving from (a) to (b), where two parcels have been swapped. (c) shows the formation of minimum energy, where the field strength and susceptibility now both decrease with height. (Online version in colour.)

    We thus look to prove the system is unstable using linear stability analysis. We consider perturbations to the stationary state such that

    where u^(z)=(u^(z),v^(z),w^(z))T, ζ=ei(αx+βy)+st, ϵ1, the prime denotes d/dz, α, β are real and positive wavenumbers, the hat variables could be complex, s is the growth rate of the disturbance and could be complex and Re denotes the real part. From now on, the real part will not be written explicitly.

    After substituting (3.7) and linearizing, (2.10) in component form is

    (2.9) and (2.11) give
    respectively, where
    We require w^=0 at z=z1,2, and it follows that χ^=0 at z=z1,2. Note that by considering the sum of (3.8) multiplied by α and (3.9) multiplied by β, we find the equations are satisfied by αv^=βu^. The direction of the wavevector (α, β) is arbitrary, and we could have chosen β=0 without loss of generality.

    ϕ^(l) has the general solution

    for constants q1(l),q2(l). Imposing ϕ(l) regular at z± gives
    Substituting (3.7) into (2.4) and linearizing gives
    Integrate (3.16) to give
    Continuity of B and ϕ at the walls, and (3.15) result in
    at z=z2, and
    at z=z1.

    Manipulating the equations, as shown in appendix A, produces the eigenvalue equation

    Equations (3.17) and (3.12) give Hϕ^=0 and (Hϕ^)=0 at z=z1,2, as well as ϕ^ satisfying (3.19) and (3.20).

    Solving (3.21) determines the eigenmodes, ϕ^, and the associated eigenvalues, s, but here we demonstrate the existence of unstable modes for general H0. Multiply (3.21) by Hϕ^, where ϕ^ is the complex conjugate of ϕ^, to obtain

    Using (3.18),
    Integration by parts on (3.23) and invoking the boundary conditions gives
    Substituting (3.24) into (3.22) and performing integration by parts on the other terms we obtain
    Equation (3.25) is of the form as2+bs+c=0, with a,b>0, c<0 and all are real. Consequently s is real. It follows that ϕ^ is real, as (3.21) is a linear equation with constant coefficients. Observe from (3.17) that χ^ is real, and therefore w^ is real, but (3.8) and (3.9) indicate that u^ and v^ are imaginary. All wavenumbers, therefore, possess a growing mode (s>0) and are unstable. This agrees with FBM that when dH/dχ<0, the system is unstable.

    4. Stabilization with a rotating field

    We seek to dampen the unstable modes. In an effort to stabilize the system, we apply a constant field in the x and y directions, such that H0=(D,E,H0(z))), but the equilibrium still requires χ0=A/H0(z)1, u0=0 and p=p0(z). Analogous analysis to §3 produces two simultaneous equations

    Observe that modes satisfying Dα+Eβ=0 remain unstable, since the equations return to that of §3. Thus, perpendicular modes to the horizontal field are not dampened and will grow.

    Applying an alternating field is analogous to applying a constant field, in that there will exist modes perpendicular to the field which remain unstable. But, for a rapidly rotating field, the field direction changes sufficiently often that an unstable mode perpendicular to the field may not have time to grow before it is no longer perpendicular to the field, and is dampened by the field. We now investigate whether a rapidly rotating field is sufficient to stabilize the system.

    Let H=(Dcos(ωt),Esin(ωt),H0(z)), and the equilibrium remains χ=1+A/H0(z), u=0 and p=p0(z) in the fluid, and H=(Dcos(ωt),Esin(ωt),A) in the plates. Performing analogous linear stability analysis to §3, but now assuming perturbations of the form ϕ^(z,t)eiαx+iβy for each variable, we obtain two simultaneous equations,

    where at z=z2,
    and at z=z1,
    By requiring u=0 at z=z1,2, we infer tχ^=0 at z=z1,2. Thus χ is constant at the walls, which we define to be zero. We set D=E for simplicity, but the analysis is analogous for DE.

    We suppose there are two time scales; one for the growth rate of the instability and the other for the rotation speed such that

    The scaling in (4.8) anticipates that χtuχ in (2.11). Boundary conditions for χ give χc(z,t)=A3,4=0 at z=z1,2. The boundary conditions for ϕ^ give
    at z=z2, and
    at z=z1.

    To obtain an eigenvalue equation for ϕ, we substitute (4.7) and (4.8) into both (4.3) and (4.4). In the resultant equations, we assume O(1/ω2) and O(1/ω) terms are negligible, equate cos(ωt) and sin(ωt) terms and time average over the fast time scale. We obtain the equations (B 5), (B 6), (B 10) and (B 12) given in appendix B alongside their derivation. Since the coefficients of said equations do not depend on t, we assume tχc=sχ(z)est, tϕc=sϕ(z)est. Manipulating (B 5), (B 6), (B 10) and (B 12) we obtain the eigenvalue equation,

    for which the derivation is shown in appendix B. Moreover, at z=z1,
    and at z=z2,

    Next, suppose D is sufficiently large to give

    such that (4.11) approximates as
    Rather than solve (4.15), we determine whether a sufficiently large D dampens all disturbances. Multiply (4.15) by
    integrate over the domain and perform integration by parts to give
    Similarly if we multiply (4.15) by
    we obtain

    Although (4.17) and (4.19) are not quadratic equations due to ϕ depending implicitly on s, we write them in the form ajsj2+bjsj+D2cj=0, where j=1,2 for (4.17) and (4.19) resepectively. s must satisfy (4.17) and (4.19) simultaneously, and observe that a1>0, b2>0 and cj>0. If a2>0 and b1>0, Re(s)<0 and the system would be stable. Using the quadratic formula and expanding sj for D1 gives

    where we use the modulus sign to denote that the variable is positive. Thus to highest order
    Since s1 is imaginary and to leading order s must equal both s1 and s2, it follows that s2 must be imaginary, and therefore a2>0. To next order
    Since a2>0, it follows that b1>0, and therefore Re(sj)<0. We conclude that for sufficiently large D, all modes are stable. We have thus demonstrated that adding a sufficiently strong rapidly rotating field in the (x,y) plane stabilizes the system.

    5. Concluding remarks

    This paper investigates the stability of an arbitrary planar equilibrium where the susceptibility and field vary in the normal direction. The configuration is proven to be unstable. Yet, a sufficiently large, rapidly rotating magnetic field will dampen all unstable modes, where a constant or alternating field cannot.

    We acknowledge that a rapidly rotating field may have thermal effects, but we assume the FF to be isothermal throughout. It should be noted that for an FF with finite electric conductivity, a time-varying field may induce an electric field, yet for an insulating ferrofluid Faraday's law remains negligible. Moreover, a rotating field may give rise to magnetoviscocity in the FF, as mentioned in §2. It is found that for weakly non-equilibrium situations (and therefore also for quasi-equilibria), a rotating field does not induce a ‘spin up’, as would an alternating field [23]. Under the quasi-equilibrium approach, the magnetization relaxation time of the particles in the FF is much less than the time scale of the hydrodynamic processes of interest [24], and therefore magnetoviscous effects are consistently ignored in this paper. For further discussion of magnetoviscous effects see for example [2,25,26] and [27].

    The overall stability of the FF configuration is of course vital for applications. For example, in a FF sealant or drug delivery mechanism, the ferrous particles can be held in place using a rotating field, and released at will by turning it off. Conveniently, here the field is generated outside the FF, whereas in FBM the equilibrium required a wire running along the axis of a FF column to produce an azimuthal field. In some instances the presence of a wire may not matter, but in others it could be a prohibitive obstruction. The results of this paper may, therefore, widen the range of potential applications of non-uniform FFs.

    More generally, when the field is a function of the susceptibility, equipotential surfaces of zero mean curvature give a set of solutions to the governing equations, and the planar geometry can be regarded as a special case of this. Thus, in theory, for a general geometry there exist stationary states where dχ/dH<0, if the boundary conditions are satisfied. FBM argues that such equilibria are unstable, since the regions of highest susceptibility and regions of strongest field do not coincide. This paper proves that for many planar configurations such instability can be suppressed by an additional rotating field. For more general three-dimensional configurations of χ, including those with zero mean curvature, a single rotating planar field may not dampen all instabilities, but two such fields should suffice to maintain stability. Establishing other equilibria and ensuring their stability could extend the possible applications of FFs still further.

    Data accessibility

    This article has no additional data.

    Authors' contributions

    S.H.F.B.: formal analysis, investigation, resources, writing—original draft, writing—review and editing; A.J.M.: formal analysis, investigation, project administration, supervision, writing—original draft, writing—review and editing.

    All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

    Conflict of interest declaration

    We declare we have no competing interest.


    This work was supported by the Engineering and Physical Sciences Research Council who awarded the Industrial Strategy Scholarship to S.H.F.B.

    Appendix A

    Apply the operator

    iα(sρη2)A 1
    to (3.11) to give
    α2(sρη2)u+αβ(sρη2)v=iα(sρη2)w.A 2
    Substitute (3.10) into (A 2) to give
    (α2+β2)(sρη2)u=iα(sρη2)w,A 3
    and take the derivative w.r.t z to give
    (α2+β2)(sρη2)u=iα(sρη2)w.A 4
    Substitute (A 4) into (3.9) to give
    (sρη2)2w^=μ0(α2+β2)(H0H0χ^+AH0H0ϕ^)A 5
    and substitute (3.12) and (3.17) to give
    (s2ρηs2)2Hϕ=μ0(α2+β2)AH0(ϕ^+H0Hϕ)H0.A 6
    Substituting (3.18) into (A 6) gives
    (s2ρηs2)2Hϕ=μ0(α2+β2)2AH00zϕ^(ζ)H0(ζ)dζ.A 7

    Appendix B

    Substituting (4.7) and (4.8) into (4.3) gives

    ωρ2[H02H0(t2χc+(ttA3ω2+2tA4ωA3)cos(ωt)+(ttA4ω2A42tA3ω)sin(ωt))]η4[H02H0(tχc+(A4+tA3ω)cos(ωt)+(A3+tA4ω)sin(ωt))]=Aμ0(α2+β2)[H0H0(χc+A3ωcos(ωt)+A4ωsin(ωt))+AH0H02(iD(αcos(ωt)+βsin(ωt))+H0z)(ϕc+A1cos(ωt)+A2sin(ωt))],B 1
    and into (4.4) gives
    (z(AH0z)AH0(α2+β2))(ϕc+A1cos(ωt)+A2sin(ωt))=iD(αcos(ωt)+βsin(ωt))(χc+A3ωcos(ωt)+A4ωsin(ωt))z(H0(χc+A3ωcos(ωt)+A4ωsin(ωt))).B 2

    Assume O(1/ω2) and O(1/ω) terms are negligible. Equating terms in (B 1) gives

    ρ2(H02A3H0)+η4(H02A4H0)=μ0A2(α2+β2)H0(iDαϕc+H0zA1)H02B 3
    ρ2(H02A4AH0)η4(H02A3AH0)=AH0μ0(α2+β2)(iDβϕ+H0zA2)H02.B 4
    Similarly, (B 2) gives
    MA1=iDαχcB 5
    MA2=iDβχc,B 6
    M=z(AH0z)AH0(α2+β2).B 7
    Time averaging over the fast time scale in (B 1) and (B 2) gives
    ρ2(H02t2χcH0)η4(H02tχcH0)=μ0A(α2+β2)(H0H0χc+AH0(2H0zϕc+iD(αA1+βA2))2H02)B 8
    Mϕc=z(H0χc).B 9
    Since the coefficients do not depend on t we assume tχc=sχ(z)est, tϕc=sϕ(z)est to give
    (ρs22ηs4)(H02χAH0)=μ0(α2+β2)(H0H0χ+AH0(2H0ϕ+iD(αA1+βA2))2H02)B 10
    Mϕ=(H0χ).B 11
    Integrating (B 11) gives
    χ=AH0HϕH02B 12
    and substituting (B 12) into (B 10) simplifies to
    (ρs22ηs4)(Hϕ)=μ0A(α2+β2)H0((α2+β2)0zϕ(ζ)H0(ζ)dζ+iD(αA1+βA2)2H02).B 13
    Apply the operator
    MH02H0B 14
    to (B 10) and substitute for MA1, MA2 and (B 12), to obtain the eigenvalue equation,
    M(H02H0(s2ρ2sη4)Hϕ)=μ0A(α2+β2)2(M(H020zϕ(ζ)H0(ζ)dζ)+AD2H0HϕH02).B 15


    Published by the Royal Society under the terms of the Creative Commons Attribution License, which permits unrestricted use, provided the original author and source are credited.