Numerical evaluation of the ferrofluid behaviour under the influence of three-dimensional non-uniform magnetic field

Ferrofluids have always been the centre of interest for a broad range of applications where the flow manipulation is achieved using different types of magnetic fields. In the present work, a Finite Volume Method based numerical study is performed to investigate the ferrofluid behaviour in the laminar flow regime using the complete set of Ferrohydrodynamics equations. Four distinct single layer finite solenoids of different L s / D s (length to diameter ratio) are considered to replicate a 3D (three-dimensional) non-uniform magnetic field, and a twofold validation is also reported to authenticate the generated field distributions. Moreover, a magnetic field dependent viscosity model is used to take into account the variation of ferrofluid viscosity in the presence of a non-uniform magnetic field. To accommodate the varying particle size distribution of ferrofluid, an effective relaxation time constant is used that can accurately predict the magnetization using Debye relaxation mechanism. The obtained flow characteristics are then further discussed and compared along with the case of stationary ferrofluid for all solenoid configurations. New insights about the flow structures are provided using the λ 2 criteria of vortex identification technique and the critical point theory analysis. It is observed that each solenoid arrangement have a distinct magnetic field distribution which leads to the occurrence of unique vortex structures and fluid deformations. Also, different flow trajectories are noticed within the fluid domain owing to the peculiar distribution of Kelvin force density. Results from the present work have direct implications in all ferrofluid based momentum transport devices and can be effectively used to further improve our understanding about their performances.


Introduction
The continuum description of FerroHydroDynamics (FHD) as first defined by Neuringer and Rosensweig (1964), deals with the interaction of fluid dynamics and thermodynamics of magnetically polarized fluids.These incompressible fluids consist of colloidal suspensions of magnetic particles with the surfactant coatings to avoid any coalescence even in the presence of considerably strong magnetic fields (Rosensweig, 1985;Rosensweig, 1987).In the absence of an external field H, despite having individual magnetic moments, the net magnetization of these particles (and eventually of the ferrofluid) is always zero, as the particles are in continuous Brownian motion.However, when subjected to a magnetic field, the particles will try to align their dipoles with the applied field according to either of the two different relaxation mechanisms (Shliomis, 1974) (Neel's and Brownian), based on the particle size, ferrofluid's viscosity, and the magneto-crystalline anisotropy constant (Odenbach and Thurm, 2002).As a result of the required relaxation time and/or the considerable flow vorticity, the magnetization of ferrofluid M always lag with respect to the applied magnetic field and this deviation leads to a body couple μ 0 M × H, known as torque density which results in the particle rotation.Moreover, in the presence of a non uniform magnetic field, magnetic particles will experience a volumetric force μ 0 (M ⋅∇) H, named as Kelvin Force Density, which drives them towards the region of higher magnetic field.Both the translational and rotational motion of these particles is then transferred to the corresponding ferrofluid by shear mechanism, spin diffusion, and rotational friction.These resulting body forces are the promising reason which make ferrofluids a preferential choice in all the applications where ability to regulate the flow is a decisive factor.This majorly includes mechanical seals, damping systems, loudspeakers (Rosensweig, 1982;Bailey, 1983), MEMS devices (Perez-Castillejos et al., 2000, 2001, 2004, 2005), and heat exchangers (Barclay, 1982;Gandomkar et al., 2017;Bezaatpour and Rostamzadeh, 2020) among several others (Rosensweig, 1985;Odenbach, 2009;Berkovski and Bashtovoy, 1996).Additionally, on the similar lines of FHD, Bio magnetic fluid dynamics is also explored by many researchers (McKay et al., 2007;Türk et al., 2014;Tzirakis et al., 2016;Imran et al., 2021), where the behaviour of blood is observed under magnetic fields for cell separation, drug targetting, and monitoring blood flows during surgeries (Haik et al., 1999;Voltairas et al., 2002;Alexiou et al., 2002;Ritter et al., 2004;Kenjereš, 2008).
Another prominent characteristics of the ferrofluids is to have a varying rheological behaviour with respect to the nature of applied magnetic field.Ferrofluids show increase in the viscosity under a steady magnetic field, and at low frequency oscillating fields, whereas they exhibit reduction in viscosity at higher frequencies (Shliomis and Morozov, 1994) (f > 1/τ B ).After the initial experimental (Moskowitz and Rosensweig, 1967) and the theoretical work (Zaitsev and Shliomis, 1969) on the spin up phenomenon, Shliomis and Morozov (1994) demonstrated that in the presence of shear flow, oscillating magnetic fields produce a non zero angular velocity, and thus reduces the overall viscosity of the fluid.This negative viscosity effect is then explored by many researchers experimentally (Zeuner et al., 1998;Schumacher et al., 2003;Shliomis and Zaks, 2006;Tekir et al., 2020), theoretically (Felderhof, 2001;Krekhov et al., 2005), and numerically (Sánchez and Rinaldi, 2010;Finlayson, 2013;Jimenez et al., 2020) to observe and control the flow behaviour.
On the similar lines, several studies have been performed in the context of steady magnetic field to examine the interplay between flow shear and pressure gradient within the magnetic fluid.The effect of externally applied transverse magnetic field on a plane Poiseuille flow is numerically compared by Gedik et al. (2012).By observing the pressure drops and flow distribution, it is summarised that flow velocity exhibits proportional reduction with respect to the applied magnetic field.Similar configuration is recently investigated by de Carvalho and Gontijo (2020) for magnetization diffusion in a parallel duct.Authors discussed about three different timescales that controls the nonequilibrium magnetization process, and concluded that it can be controlled with the help of flow rate, field strength, volume fraction, and particle size distribution.A two phase circular Poiseuille flow is analytically studied by Recebli and Kurt (2008) and it is demonstrated that the perpendicular magnetic field, leads to the reduced local velocity for both the phases.Hart (2002) investigated a rotating Couette flow in a rectilinear geometry for its effect on magnetic viscosity.A radial field gradient is imposed using an array of permanent magnets on a vertical long channel filled with the ferrofluid.Using perturbation analysis, it is revealed that the elemental roll instability of the rotating Couette flow is suppressed by up to 10 % wavy flow behaviour.A Couette-Poiseuille flow in 2D parallel plates with a perpendicular magnetic field is numerically examined by Ghosh and Das (2019).Authors discussed the possibility to control the flow separation with the help of an adjustable magnetic field and provided expressions to predict the spin field, velocity, and wall friction coefficient.Yang and Liu (2020) studied the same configuration numerically for the influence of magnetic relaxation using more accurate magnetization equation.It is shown that the application of magnetic field decreases the volume flow rate and it further reduces with the increase in the field strength before achieving its saturation value.
Moreover, because of its applicability and relevance in several applications (Bahiraei and Hangi, 2015;Kole and Khandekar, 2021), many researchers extensively analysed ferrofluid flow in 3D ducts/pipes of circular cross-sections.Krekhov and Shliomis (2017) numerically and analytically investigated the core rotation of ferrofluid flow in the presence of an uniform axial magnetic field.Authors proposed the existence of two different phases where the inner core follows swirling path and outer annulus shows pure axial characteristics.Vafeas et al. (2019) numerically examined the effect of a transverse field on a 3D curved circular duct using continuity vorticity pressure variational method.It is mentioned that, the axial pressure drop varies approximately in linear manner with the applied strength curvature and an external field result in the formation of secondary flows.Papadopoulos et al. (2012) carried out an extensive numerical work of examining the effect of several direct current solenoid configurations on a laminar ferrofluid flow.Authors concluded that, the steady magnetic field alters the axial pressure distribution inside the pipe and this pressure drop varies linearly with the volume fraction.Some similar attributes are also observed for ferrofluid duct/pipe flows under the influence of time varying magnetic fields.An experimental and numerical work carried out by Schumacher et al. (2003) proposed that flow rate, field strength, and frequency of the magnetic field affects the measured pressure drop within a pipe flow.Based on the validation with the aforementioned experimental results, Krekhov et al. (2005)  studied numerically by Felderhof (2001) for comparing the influence of three different magnetization equations, and concluded that there is severe discrepancy among them at higher amplitudes.Mao et al. (2011) studied a closed loop pumping mechanism in the absence of any moving components using the spatially travelling magnetic fields.On the basis of experimental and the numerical observations, authors mentioned about the possibility of reversible dimer formations at moderate field magnitudes and subsequent pumping effect even with small volume fractions.
As mentioned previously, there have been several studies which depicts the application of steady magnetic field (Krekhov and Shliomis, 2017;Papadopoulos et al., 2012;Vafeas et al., 2019;Asfer et al., 2016;Kaneda and Suga, 2018) or time varying magnetic field (Krekhov et al., 2005;Felderhof, 2001;Schumacher et al., 2003;Mao et al., 2011;Tekir et al., 2020;Jimenez et al., 2020) in the circular pipes of straight (Krekhov et al., 2005;Krekhov and Shliomis, 2017;Papadopoulos et al., 2012;Felderhof, 2001;Schumacher et al., 2003;Asfer et al., 2016;Kaneda and Suga, 2018;Tekir et al., 2020;Jimenez et al., 2020) or curved cross-sections (Vafeas et al., 2019;Mao et al., 2011).Nevertheless, there is no work focussing on the flow behaviour of ferrofluids in low Re regime, with 3D non-uniform magnetic fields.Such small Re flows are often experienced in several applications such as ferrofluid actuators (Yamahata et al., 2005), micro-pumps (Hatch et al., 2001(Hatch et al., , 2004)), and MEMS devices (Perez-Castillejos et al., 2000).Thus, present work aims to examine ferrofluid characteristics at three different magnitudes of Re (1, 25, and 50), along with stationary case (Re = 0), each under the influence of four non uniform field distributions generated from the single layer finite solenoid.Additionally, although Papadopoulos et al. (2012) done a pivotal work of investigating ferrofluid behaviour with finite solenoids, no details are mentioned about the variation and the influence of magnetic field distribution in radial direction.Thus, in the present study, both radial and axial field distributions are discussed in depth for the considered four different solenoid configurations.The numerical results obtained from the present validated model has direct implications in all the fields where ferrofluid is used for the momentum transport.
The present work is organised as follows, Section 2 includes information about computational domain, governing equations, boundary conditions, and thermophysical properties, whereas in Section 3 detailed numerical methodology (field distribution of solenoids, grid independence study, model validation, and model specification) is mentioned.After discussing about results in Section 4, a concluding remark is given in the Section 5.

Geometrical details
In the present work, a long horizontal pipe of circular cross-section of a hydraulic diameter D h and length L p equivalent to 15 D h is considered, which is completely filled with the ferrofluid.The present geometry is specifically chosen as it is practically relevant and is widely used in several ferrofluid applications such as pumping devices and actuators.The four different solenoid configurations which are considered in the present study are mentioned in ascending L s /D s in the Table 1, based on their distinct length and diameter.These single layer finite solenoids are positioned coaxially with the pipe, and are used to produce the steady, non-uniform magnetic field distributions as mentioned in 3.1.A Cartesian coordinate system (x, y, z) is used for most of the work (and explicitly mentioned if otherwise), where the axial direction is defined along the z axis as shown in Fig. 1.

Governing equations and boundary conditions
To determine the time dependent behaviour of a laminar, incompressible, Newtonian, and an isothermal ferrofluid, standard equations of FHD are considered.The Maxwell's equation, continuity equation, and the momentum equations (both linear and angular) are solved along with the magnetic relaxation mechanism.Owing to their insignificant electrical conductance (Neuringer and Rosensweig, 1964;Rosensweig, 1985;Shliomis, 2002), the Maxwell's equations are solved in its magneto-statics formulation (Jackson, 1999;Rosensweig, 2004) as shown in Eq. 1.
Maxwell's Equations ∇⋅B = 0, ∇ × H = 0, and B = μ 0 (H + M) (1) Here, B is magnetic induction vector, H is magnetic field vector, and both are coupled via magnetic permeability (μ 0 ) and magnetization vector M.For examining the flow characteristics and its interaction with the externally applied magnetic field, the continuity equation is considered along with the linear momentum equation and angular momentum equation as shown in Eqs. ( 2)-( 4) respectively (Dahler and Scriven, 1961).Continuity Equation As discussed earlier, the magnetic particles can have distinct behaviour under the influence of an externally imposed magnetic field.Thus, to predict such unique responses of the magnetic particles and eventually of the ferrofluids, several magnetization equations have been proposed (Shliomis, 1971;Martsenyuk et al., 1974;Shliomis, 2001;Rosensweig, 1985;Felderhof and Kroh, 1999;Ivanov and Kuznetsova, 2001).Based on its accuracy, ease of computation, and range of validity (Ωτ eff < 1), a phenomenological magnetization equation (Shliomis, 1971) is used for the present work, which is generalisation of the Debye relaxation (Debye, 1929) as shown in Eq. 5.
Magnetization Equation

∂M ∂t
In Eqs. ( 2)-( 5), U represents the velocity vector, S is the internal angular momentum vector, and p is the pressure field.Also, τ eff , and τ S are used to define the effective relaxation time of particles and spin relaxation time, where I is the moment of inertia of particles per unit mass.Additionally, ρ nf denotes the density of the colloidal suspension, whereas μ eff and η ′ refer to the dynamic viscosity and the spin viscosity offered by the ferrofluid.Moreover, M 0 is the equilibrium magnetization which is defined using the Langevin function L(ξ) as shown below In Eq. 6, n is the number density of particles (n = 6ϕ/πd 3 ) defined for a specific particle diameter (d) and volume fraction (ϕ), whereas K B is the Boltzmann constant, and H corresponds to the mean magnetic field ) experienced by the ferrofluid.Also, m is the magnetic moment of a single particle of Fe 3 O 4 which have a face-centered cubic structure with a coordination number of 8, where each molecule possess a magnetic moment of 4μ B (Kittel et al., 1996), and thus m is given as (Aminfar et al., 2013) Now, for tiny magnetic particles (of diameter ≈ 10 nm), the magnitude of I is insignificant ( of the order 10 − 17 m 2 ) (Schumacher et al., 2008), hence the total derivative in Eq. 4 can be neglected (Shliomis, 2001;Schumacher et al., 2003;Finlayson, 2013;Ghosh and Das, 2019).
Furthermore, owing to the very small value of η ′ (10 − 17 − 10 − 18 ) (Feng et al., 2006;Rosensweig, 2013), contribution from the diffusion term is also insignificant.As a result of this, Eq. 4 simplifies to After substituting from above relation, the linear momentum Eq. 3 can be modified as The last two terms of Eq. 9 represent the body forces experienced by the ferrofluid, known as Kelvin Body Force (also known as Kelvin Force Density (KFD)) and torque density (also known as Magnetic Torque (MT)) respectively.As mentioned earlier, KFD arises in the ferrofluid because of the non uniform distribution of H, whereas MT refers to the torque experienced as a result of the difference between flow vorticity and internal spin of ferromagnetic particles.Similarly, after replacing S (from Eq. 8) in the magnetization Eq. ( 5) and then using the relations for τ S (ρ s d 2 /60μ eff ),I (ρ s ϕd 2 /10), and ϕ (nπd 3 /6), the Eq. 5 can be rewritten as ∂M ∂t In general, because of the non-uniform particle size distribution within a ferrofluid (Kaiser and Rosensweig, 1969;Rosensweig, 2013), it is possible that relaxation can occur simultaneously by both (Neel's and Brownian) mechanisms (Mao et al., 2011).Thus, for the present work an effective magnetic relaxation time constant (Martsenyuk et al., 1974) is used as mentioned below Here, f 0 and K represent the Larmor frequency and the magnetocrystalline anisotropy constant respectively (Table 2), whereas V is the core volume of the solid particle alone ) , and V h refers to the overall hydrodynamic volume of the nanoparticle including surfactant coating ) .Thus, for every combination of H and T, we have a specific magnitude of M 0 and the ferrofluid will try to approach this corresponding equilibrium value based on Ω, H, ϕ, μ eff , and τ eff .These aforementioned set of governing equations of the present study are then put to a closure along with the magnetic boundary conditions (refer Appendix A) and solved numerically using the initial and the boundary conditions mentioned in Eq. 12, where ρ is defined as a radial coordinate , and n refers to the normal direction.

Non-dimensionalization and thermophysical properties
To reduce the number of variables and make the present analysis generically applicable, current governing equations are scaled with the a See Supplementary Material of Mao et al. (2011).
help of dimensionless parameters mentioned in Eq. 13 Here, D h is the hydraulic diameter of the pipe, and U refers to the mean velocity magnitude.Additionally, M 0,max is the maximum equilibrium magnetization corresponding to ϕ = 1, and a sufficiently large value of magnetic field (H = H max ) i.e., L(ξ) = 1.Also, H max is the maximum strength of magnetic field (Eq.17) obtained from an electric wire having N number of turns over solenoid length L s , and carrying a steady current of fixed magnitude I through it.After substituting the corresponding values from Eq. 13, the governing equations (Eq.2-5) can be written in the dimensionless form as shown below, where, Similarly, the initial and the boundary conditions are also nondimensionalised and are written in the modified form as shown below, in which both ρ and z are scaled with the hydraulic diameter D h .
In the present study, a colloidal suspension of water based Fe 3 O 4 nanoparticles is considered as a working fluid.Thus, using corresponding ϕ, individual contribution of the base fluid and the solid particles (Table 3) is considered while calculating the thermophysical properties.The density (ρ nf ) for the present ferrofluid are calculated as per the Eq.19 as shown below Additionally, for describing the flow resistance in the presence of a magnetic field, an effective viscosity model is used which is proposed by Shliomis (1971).As expressed in Eqs. ( 20)-( 22), it takes into account the enhancement of viscosity experienced by the ferrofluid in response to the applied magnetic field.
In Eq. 20, μ nf is calculated with the help of a specific correlation proposed by Rosensweig (1985).It provides additional corrections to the well established Einstein's viscosity model (Einstein, 1905) to consider the behaviour of ferrofluids at higher ϕ.
Here, ϕ c is the critical value of concentration beyond which the results obtained from Eq. 21 differs with the experimentally measured values.
For the present case, its value is taken as 0.74, since it provides accurate viscosity value at all higher volume fractions up to ϕ = 30% (Buschow, 2003).The additional increase of viscosity in the presence of an external magnetic field, often termed as rotational viscosity is given as In Eq. 22, β is the angle between two vectors Ω and H, whereas d and s refers to the particle diameter and surfactant thickness, which for the present study are considered as 10 nm and 2 nm respectively (Buschow,

Table 3
Thermo-physical properties of base fluid and nanoparticles used (Sun et al., 2020).
S. Dalvi et al.

2003)
. Moreover, ξ denotes the ratio of the magnetic energy and the thermal energy of the particles as described in the latter part of the Eq. 6.

Magnetic field of a finite solenoid
One of the foremost work on the magnetic field distribution of thin, finite, and a circular cross-sectional solenoid is carried out by Callaghan and Maslen (1960) using magnetic vector potential.Authors derived radial and axial field variations for different lengths and diameters of solenoid using complete elliptic integrals.Several other studies (Conway, 2001;Dasgupta, 1984;Derby and Olbert, 2010;Muniz et al., 2015;Labinac et al., 2006) explain these field distributions in distinct manner using different methods and analytical equations.One of many such works is done by Labinac et al. (2006), in which authors discuss about theoretical expressions for the field distribution of a single layer solenoid using various mathematical functions such as the Bessel function, hypergeometric function, and the Legendre polynomials.The aforementioned theoretical expressions are written in the cylindrical coordinates for an azimuthally symmetric solenoid which axially starts from the origin.These field distributions are then generalised by Papadopoulos et al. (2012) for a thin solenoid of finite length which can be located arbitrarily anywhere in the axial direction.The generalised distribution is defined using the Bessel functions of the first kind for a solenoid diameter D s and solenoid length L s (which extends from z 1 to z 2 such that z 2 − z 1 = L s ) with N number of wire turns that carries electric current of magnitude I as shown in Eq. 23-25 where f(k; z) is defined as a piece-wise function based on the value of z as shown below The radial and the axial field distribution obtained from the above . For all four configurations, both H ρ and H z have maximum magnitudes in the vicinity of the solenoid (at 2ρ/ D s = 1).Additionally, it can be clearly seen from Fig. 2, at the solenoid surface and near both the edges of the solenoid, the magnitude of H ρ tends to infinite which reflects from the sharp peaks at Ls start and Ls end .Moreover, at smaller L s /D s , there exists a discrepancy between the magnitude of H ρ at both the ends (sub-Fig.2a and b where , and such a peculiar behaviour can also be confirmed from the results of Callaghan and Maslen (1960).Further, it is depicted in Fig. 3 that, for a specific solenoid length, increase in D s leads to the reduced magnitudes of H z , whereas at a particular solenoid diameter, increase in L s leads to ascending H z .It can also be noticed that the radial variation of H z (from 2ρ/D s = 0 to 2ρ/D s = 1) decreases with increase in L s /D s .

Grid distribution and grid independence study
A non-uniform and a non-orthogonal, structured grid is implemented for the present analysis, in which the local grid refinement is done with a linear expansion rate.Results from the conventional orthogonal grid distribution (O-grid) reveals the presence of artificial flow disturbances at the diagonal locations in the streamwise direction.This supplementary discretization error is attributed to the block structured arrangement of orthogonal grids for the circular domains, and thus avoided by implementing a non-orthogonal body fitted grid.Additionally, as a consequence of the no slip condition, fine resolution is provided near the solid wall as shown in Fig. 4a to consider the boundary layer effects.As it can be clearly seen from Fig. 2 and Fig. 3 that for all four solenoid configurations, the field magnitude have large gradients near both the ends of the solenoid.Thus, to capture their influence accurately, grid elements with smaller spacing are also specified near both the ends in the axial direction as shown in Fig. 4c.
Moreover, to ensure the least possible discretization error, 7 different grid configurations with distinct spacing arrangements are compared for their influence on the flow variables as shown in Table 4.The assessment is done for the S4 solenoid configuration with Re = 1 and Rm = 1 × 10 8 .
The surface averaged mean velocity magnitude ( Ũ) as shown in Eq. 26 is compared at a radial mid section (z = 7.5D h ) and an axial mid section (x = 0).
Based on the obtained results, Grid 2 is finalised for all the remaining computations of the present study and its distribution is shown along two cross-sections in Fig. 4. For this finalised grid, radial and the axial spacing between two consecutive grid elements (Δ N) is also shown in Fig. 4b and Fig. 4d respectively.However, it is to be noted that, irrespective of the different spacing arrangement, the total number of grid elements are kept same for both L s = 1 D h and L s = 5 D h .

Model validation
To verify the accuracy of the magnetic field distributions (from the Section 3.1), the present numerical model is validated with both the experimental measurements and the theoretical work of Labinac et al. (2006).The comparison is carried out for a specific configuration of single layer solenoid of length 279 mm and radius 150 mm, having 150 turns of copper wire which carries a steady current of 9 A. The radial and axial magnitudes of magnetic field are compared at a specific axial locations, as mentioned in the respective sub-figures of Fig. 5.In the context of present study, the solenoid is wrapped around a pipe and we are mainly concerned about the field distribution inside the solenoid.Thus, the comparisons are only done for the radial locations which lies within the periphery of solenoid (up to 2ρ/D s = 1).Additionally, to further compare the precision of the Bessel function calculations in the OpenFOAM framework, the field distributions are also compared against the results obtained with the help of complete elliptic integral formulations (provided by Callaghan and Maslen (1960)) using MATLAB (refer Appendix B).Lastly, it can be confirmed from Fig. 5 that, the numerical results using the Bessel functions accurately follows the work of Labinac et al. (2006) as well as of the Callaghan and Maslen (1960).
Further, to assure the correct numerical implementation of the magnetization equation (Eq.10), the present model is also validated with the numerical work of Papadopoulos et al., 2012.The nondimensional axial pressure gradient within a long circular pipe (20D h ) is observed for a laminar ferrofluid at various solenoid configurations.

The validations is done at
, ξ max = 2 × 10 − 2 , and ϕ = 0.1, where ξ max corresponds to H max as described in Eq. 17.The geometrical details, grid distribution, discretization schemes, and the solution algorithm are all kept in accordance with Papadopoulos et al. (2012).As it can be seen from the Fig. 6, results from the present numerical work agrees perfectly with those from Papadopoulos et al. (2012).Thus, it can be evidently said that, the current numerical model can accurately follow the behaviour of ferrofluid in the presence of a 3D non uniform field distribution generated from a finite solenoid.

Numerical method
The present set of governing equations are solved with the help of a C++ based open-source CFD Toolbox OpenFOAM (Weller et al., 1998), which use the Finite Volume Method based discretization for all the flow variables.The convective terms of all the governing equation are discretized using upwind biased second order accurate scheme, whereas the diffusive terms are solved using second order central difference scheme.In addition, the temporal marching of the solution is carried out with the help of second order backward difference scheme.The pressure velocity coupling is achieved using the PIMPLE algorithm (Greenshields, 2015), which is a combination of PISO algorithm and SIMPLE algorithm (Ferziger et al., 2002;Moukalled et al., 2016).It is mainly used for computing the complex transient flows, as it allows to have the larger time steps without destabilizing the solution (Holzmann, 2016).Further, while solving iterative matrices, the diagonal dominance is achieved by linearising the source terms of magnetization equation (Eq.10) using the Picard's method (Ferziger et al., 2002;Patankar, 2018).Finally, the symmetric matrices of the resulting algebraic equations are solved with the help of generalised algebraic multi grid solver, whereas smooth-Solver is used for asymmetric part along with diagonal based incomplete LU preconditioner.

Results and discussion
In the present work, the numerical computations are carried out for the four different L s /D s having H max corresponding to Rm = 1 × 10 8 and at three low Re magnitudes (Re = 1, Re = 25, and Re = 50).To further understand the sole effect of inertia forces, an additional case study of stationary ferrofluid is also observed at all L s /D s .Overall, it is observed that the 3D magnetic field distribution leads to the presence of different body forces (namely KFD and MT) which majorly affects the ferrofluid characteristics.Although, it is noticed that the magnitude of magnetic torque density is 10 − 7 − 10 − 9 times smaller than that of KFD, and thus have a negligible effect.This finding is in alignment with the results reported by Papadopoulos et al. (2012), where it is mentioned that the torques density is only influential in the case of uniform magnetic field.Thus, in the subsequent sections, only effect of KFD is mentioned while discussing the flow behaviour.

Bulk flow topology
For the current study, the global flow structures are observed by visualizing the regions of vortex formations, as vortices are majorly associated with the significant variations inside the fluid domain.At lower Re however, the vortex cores can not be identified using the conventional pressure minimum principle (Jeong and Hussain, 1995).Thus, for the present study, λ 2 criteria of vortex identification technique is implemented in which the second invariant of the local pressure Hessian (p ij ) is used.
After neglecting the unsteady irrotational straining and the viscous effects (as they lead to the insufficient criteria of pressure minima at the vortex cores) from Eq. 27, the local pressure extremes are evaluated with the help of negative eigenvalues of symmetric tensor S ik S kj +W ik W kj which is shown below ) and The symmetric part S ij and the anti-symmetric part W ij of the velocity gradient tensor (A ij ) represents the rate of strain and the rate of rotation respectively.With the help of these tensors, the vortex structures are visualized and are compared for the iso-surfaces corresponding to the value of λ 2 = − 0.01 at stationary Re (for different L s /D s ) and at smallest L s /D s (for different Re) as shown in Fig. 7 and Fig. 8 respectively.In general, it can be clearly noticed from both the figures that the major topological changes occur in the vicinity of the solenoid ends.This is completely attributed to the presence of large magnetic field gradients near Ls start and Ls end , which results in the existence of significant KFD around the solenoid ends.Further, owing to the absence of any angular variation in the present problem description, the observed vortices exhibit azimuthal symmetry for all the configurations.Although both radial and axial magnetic field gradients are responsible for the obtained vortex formations, but it is noticed that H z have a significant influence as compared to H ρ .
For the dissimilar solenoid arrangements (L s /D s ), distinct vortex structures are observed at different locations for Re = 0 as shown in Fig. 7.In this case, the lack of inertia forces results in the development of axisymmetric vortex distributions (relative to z = L s /2).For Re = 0, KFD is the only influential force experienced by the ferrofluid, thus to clearly interpret the effect of L s /D s , the obtained vortices are represented in terms of non-dimensional axial KFD component (KFD + z ).The scaling   and S3), the singular "toroidal" vortex ring is observed at Ls start and Ls end , whereas "thin cylindrical ring" is obtained for closely spaced solenoid arrangements (S2 and S4).The obtained structures can be explained as, for S1 and S3, the ferrofluid is subjected to a much wider gradient of H z in the axial direction as compared to the S2 and S4 arrangement.This can be confirmed from Fig. 3, where H z have significant non-zero magnitude across 0.2 < Z * < 0.8 for S1 and 0 < Z * < 1 for S3.Such broad magnetic field distributions result in the existence of KFD gradients over vast length, hence creating an expanded toroidal vortex at both the solenoid ends for λ 2 = − 0.01.On the other hand, for S2 and S4, H z exhibit reasonable magnitude only within 0.4 < Z * < 0.6 and 0.2 < Z * < 0.8 respectively.This narrow distribution restricts the domain of influence for corresponding KFD, thus resulting in the formation of thin cylindrical rings.Additionally, the distribution of KFD + z for λ 2 = − 0.01 have a positive magnitude at Ls start and an equivalent negative value is observed at Ls end for all solenoid configurations.The magnitude of KFD + z at both the solenoid ends also demonstrates an increment with increase in length to diameter ratio of solenoid.These KFD + z characteristics are attributed to the axisymmetric distribution of H z , as can be confirmed from Fig. 3 which also shows that the magnitude of H z is proportional to the L s /D s .
Furthermore, to analyse the effect of inertia forces, the iso-surfaces of λ 2 at different Re are also compared for the smallest solenoid configuration S1.For the case of fixed L s /D s , where KFD have identical magnitude at all Re, the obtained vortices are represented in terms of non-dimensionalised streamwise vorticity.It can be noticed from Fig. 8 that, for all Re, distinct variations of the toroidal vortex are observed at Ls start and Ls end .The obtained vortex formations are the combined outcome of both inertial forces and KFD forces that are acting on the ferrofluid.At smaller Re (Re = 0 and Re = 1), KFD is significantly dominant over its counterpart, hence axially symmetric toroidal vortex rings are observed.However, with the increase in inertia forces, KFD can no longer sustain the axisymmetric distribution of the vortices.The influence of increasing inertia force can be clearly observed from the proportional extension of vortex rings in the streamwise direction for Re = 25 and Re = 50.Additionally, as a consequence of higher inertial forces, no vortex structure is observed along the axial direction over z < Ls start for Re = 25 and Re = 50.Moreover, for all four considered cases, the magnitude of streamwise vorticity grows with increase in Re and have comparatively larger values along the central core.This augmentation around the central core is justified as the velocity distribution have the largest magnitude along the axial centreline, which leads to the occurrence of enhanced fluid movement.

Local flow characteristics
For accurate interpretation of the complex flow field, it is mandatory to have an association between the local and the global flow characteristics.Thus, to obtain it, the critical point theory is used in the present work, which is based on the collective work of Perry and Chong (1987) and Chong et al. (1990).In this method, the three invariants (P A ,Q A , and R A ) of the rate of deformation tensor (A ij ) are used to distinguish the flow topology within the fluid domain.These invariants entirely consist of the S ij and the W ij tensor (refer Eq. 28) as shown below For the present context of incompressible ferrofluid, remaining two invariants (Q A and R A ) are explored for flow interpretation in different cases.In the above formulation, flow dissipation is represented by Q A , whereas R A refers to the influence of the source/sink terms.Hence, the Joint Probability Density Function (PDF) of these two invariants are employed to decipher the complex flow field patterns.The comparisons are carried out to evaluate the influence of different L s /D s at the stationary ferrofluid (Re = 0) as shown in Fig. 11.Additionally, to acknowledge the influence of each term at all local sections, the sampling is performed specifically at three distinct radial iso -surfaces corresponding to the locations 2ρ/D h = 0.1, 2ρ/D h = 0.5, and 2ρ/D h = 0.9.Also, for improved comprehension, the aspect ratio of all the subfigures is kept same, and scaling of the invariants is done as shown below To have a clear distinction among the different flow topologies, area under the Q A − R A plot is majorly divided in the four categories as shown in Fig. 9.They are primarily described as Unstable Focus Stretching  It can be clearly observed from Fig. 11 that, for all L s /D s , the magnitude of Q * A and R * A reduces with the increase in 2ρ/D h .This distinct behaviour signifies the increase of flow rotation (Q w ) as we move towards the pipe surface at any specific solenoid arrangement.Further, it is noticed that the distantly spaced configurations (S1 and S3) exhibit an axisymmetric distribution around R * A = 0 at all 2ρ/D h .Both S1 and S3 indicate a combination of SN/S/S flow topology and UN/S/S flow topology at 2ρ/D h = 0.1 and 2ρ/D h = 0.9, whereas SFS and UFS arrangements are observed at 2ρ/D h = 0.5.On the other hand, closely spaced configurations (S2 and S4) demonstrates a SFS dominated flow structure at 2ρ/D h = 0.5.In addition to this, a clear preference for SN/S/ S structure and UN/S/S structure is observed at 2ρ/D h = 0.1 and 2ρ/ D h = 0.9 respectively.After comparing the relative magnitudes of Q * A and R * A at all radial iso-surfaces, it can be evidently put forward that, the magnetic source terms have a significant influence on the flow topology as compared to the dissipation effect.
Furthermore, the combined influence of the inertia forces and the magnetic body forces can lead to the occurrence of different deformations within the ferrofluid.Thus, to understand the specific nature of the obtained fluid deformations, the second invariant of A ij is also investigated.As it can be confirmed from the Eq.29 that, Q A represents the contributions from the S ij and W ij tensors.Thus, these tensors are used to find the relative dominance of the strain rate deformations and the rotation rate deformations by plotting the Joint PDF's of S ij S ji and − W ij W ji .As shown in Fig. 10, the distribution near to the vertical axis represent the pure strain motion, whereas the solid body rotation is observed along the horizontal axis (Soria et al., 1994).The comparison of these fluid deformations is done at a specific solenoid arrangement S4 for the different Re.It can be clearly noticed from Fig. 12 that for all Re, the magnitude of S ij S ji and − W ij W ji increases with the increase in 2ρ/D h .This reflects the rising strength of fluid deformations of each type as we shift towards the solenoid surface.Similar behaviour of growth in ferrofluid movement can also be confirmed from Fig. 11 which notify about the increase of flow rotation at higher radial iso-surfaces.At higher inertial forces (Re = 25 and Re = 50), all radial iso-surface indicate the existence of both pure strain motion and solid body rotation.Contrary to this, at smaller Re (Re = 0 and Re = 1), pure strain motion is observed only along the central core of the ferrofluid 2ρ/D h = 0.1, which is later accompanied with solid body rotation at higher 2ρ/D h .In the end, it is also noticed that the magnitude of fluid deformations increase in proportion with Re.

Distribution of volumetric forces
As it is being established for the present case that the magnetic body ] as shown in Fig. 13.Additionally, at each cross-section, every KFD component is calculated at three different radial lines (2ρ/D h = 0.1,2ρ/ D h = 0.5, and 2ρ/D h = 0.9) to capture all the variations ranging from the central core till the pipe surface.These individual components from different radial locations are then plotted in a same polar grid, to have a superior sense about the comparative influence of each component at that particular cross-section.In addition to this, the obtained magnitudes are scaled (and are represented by the superscript +) with the largest possible KFD corresponding to the maximum equilibrium magnetization ( ⃒ ⃒ M ⃒ ⃒ = M 0,max ) and maximum field gradient (|∇H| = H max ).A similar procedure is also followed to investigate the effect of Re at a particular L s /D s .However, no significant variation is observed at aforementioned radial lines, and thus for the sake of brevity these polar plots not discussed here.
In general, for all four L s /D s and at all three x − y planes, the KFD + components at 2ρ/D h = 0.9 are showing the maximum magnitude fol-Fig.13.Variation of KFD + at Re = 0, Rm = 1 × 10 8 , and ϕ = 0.05 for solenoid configuration S1 (a-c), S2 (d-f), S3 (g-i), and S4 (j-l).The respective line colours rose, green, and blue represents the x, y, and z components of KFD + .The line patterns dash dot dot, dashed, and solid in each sub-figure corresponds to the radial locations 2ρ/D h = 0.1,2ρ/D h = 0.5, and 2ρ/D h = 0.9 respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)lowed by 2ρ/D h = 0.5 and 2ρ/D h = 0.1.This observation is in alignment with the magnetic field distributions explained in Section 3.1 where higher values are noticed as we move near the solenoid surface (refer Fig. 2 and Fig. 3).Also, for closely spaced configurations (S2 and S4), every KFD + component have higher magnitudes at each cross-section than their corresponding counterparts (S1 and S3 respectively).Moreover, due to the positive gradient of H z near z = Ls start , KFD + z for all L s / D s have a positive "Omnidirectional" polar pattern at all three radial locations (2ρ/D h = 0.1, 2ρ/D h = 0.5, and 2ρ/D h = 0.9) and identical shapes with negative sign are observed for z = Ls end .On the similar lines, the inherent absence of H z gradient (refer Z * = 0.5 in Fig. 3) leads to the negligible values of KFD + z at z = Ls/2 for all solenoid arrangements.
Contrary to this, for radial components, among all three x − y planes, KFD + x and KFD + y at z = Ls/2 displays larger magnitude for all solenoid arrangements except S4.This is explained by the fact that at z = Ls/2, H ρ exhibit a point of inflection (refer Z * = 0.5 in Fig. 2) for the first three L s /D s .Thus, a significant change of magnetic field gradient exist around it, which eventually leads to the higher magnitudes of radial KFD + components.However, as L s /D s increases, this inflection point tends to becomes a stationary point over a straight line, and such an absence of gradient results in the smaller magnitudes of KFD + x and KFD + y at S4. Additionally, at all x − y planes, both KFD + x and KFD + y exhibit the characteristic "Cardioid" patterns for 2ρ/D h = 0.9, whereas "Subcardioid" and "Omnidirectional" polar shapes are observed for radial locations 2ρ/D h = 0.5, and 2ρ/D h = 0.1 respectively (in few sub-figures [namely (c), (f), (g), (i), and (l)], these aforementioned patterns are difficult to identify at first glance, as range of the radial axis is kept larger to accommodate KFD + z ).To be further specific, the "Cardioid" patterns (at 2ρ/D h = 0.9) are vertically aligned for KFD + x and are horizontally aligned for KFD + y representing the higher magnitudes in corresponding directions.However, at the smallest radial location (2ρ/D h = 0.1), for all four L s /D s , the magnitude of H ρ ≈ 0 (refer Fig. 2), resulting in the "Omnidirectional" polar patterns with negligible magnitude.And the transition between these two extreme radial locations is reflected from the "Subcardioid" shapes at 2ρ/D h = 0.5 for both KFD + x and KFD + y .

Flow field variations
To further understand the ferrofluid behaviour, the obtained flow trajectories are observed for different Re and L s /D s by visualising the streamline distributions.The individual influence of the inertia forces and the solenoid configurations are examined by plotting the flow paths for distinct L s /D s at Re = 0 and for different Re at S4.The streamlines are evaluated at x mid-plane and y mid-plane.However, owing to the angular symmetry of the present case, both planes are found to have identical flow paths.Thus, to keep the present work concise, only streamlines at x mid-plane are reported here.Moreover, it is also noticed that, for all non zero Re, the inertial forces surpass the KFD at x midplane, which leads to the conventional laminar like streamline distribution of ferrofluid at S4 configuration.For accurate comparison among all cases, the streamlines are evaluated at the same set of coordinate points which are equally spaced along axial direction for x mid-plane.To further support our understanding of the flow trajectories, the alignment of H is also visualised at x mid-plane for all solenoid arrangements.As seen from Fig. 14, the iso-magnetic field lines completely follow the divergence free behaviour in accordance with the Gauss's law (Eq.1).Additionally, for closely spaced L s /D s (S2 and S4), comparatively sharper curvatures are observed around Ls start and Ls end indicating the presence of higher field gradients.In addition to this, the proportional increase of H * z magnitude with L s /D s is also in agreement with the axial field distribution discussed in Section 3.1.
In general, it can be noticed from Fig. 15 that KFD components have the ability to notably modify the ferrofluid characteristics.This is reflected from the presence of several local recirculation zones at multiple positions for each L s /D s .Owing to the absence of any other volumetric force except KFD, the streamlines exhibit axisymmetric distribution at Re = 0 for all solenoid configurations.In addition to this, major variations in the streamlines are observed around both the solenoid ends, and are attributed to the presence of large magnetic field gradients near Ls start and Ls end .It is also observed that closely spaced L s /D s have quantitatively more recirculation zones as compared to the distantly spaced solenoid arrangements.This is mainly because of the presence of sharper curvature gradients for S2 and S4 as shown in Fig. 14b and d, which results in the likewise steep distributions of KFD across x midplane.Such KFD arrangements have a profound influence on the ferrofluid behaviour, and results in the occurrence of added recirculation zones as compared to its counterpart solenoid arrangements.Similar to this, sparsely distributed streamline patterns are noticed for S1 and S3, whereas much denser streamline patterns are observed for closely spaced L s /D s which indicates the higher magnitudes of local flow velocity.
To further evaluate the flow field variations, the distribution of normalised pressure gradient (∂p * /∂n * ) is also observed in each direction for all L s /D s at Re = 0 as shown in Figure 16.The comparison of ∂p * /∂x *  f), (i), and (l)] is performed at the central radial iso-surface corresponding to 2ρ/D h = 0.5.Overall, it is found that the magnitude of ∂p * /∂n * in all three directions is higher for closely spaced arrangements (S2 and S4) than their counterpart configurations (S1 and S3).A similar observation is previously reported by Papadopoulos et al. (2012) in the context of axial pressure gradient and for the case of non zero Re.Moreover, it can be clearly seen from Fig. 16 that at all L s /D s , both ∂p * /∂x * and ∂p * /∂y * display exactly equivalent magnitudes with a phase shift of 0.5π.This behaviour can be affiliated with the similar angular discrepancy between KFD + x and KFD + y from Fig. 13.Also, at any particular L s /D s , comparatively larger magnitudes are observed for axial pressure gradient, which depicts the significant influence of H * z on the ferrofluid behaviour.In addition to this, it is noticed that ∂p * /∂z * exhibit an increment in the magnitude at z = Ls start followed by an equivalent reduction at z = Ls end .This distinct behaviour is attributed to the positive and the negative gradient of H z at the subsequent solenoid ends (refer Fig. 3) and is in agreement with the experimental measurements of Kamiyama (1996) and the numerical results of Papadopoulos et al. (2012).

Concluding remarks
To enhance our understanding about the functioning of micropumps, actuators, and the MEMS devices, the low Re behaviour of ferrofluid is numerically investigated in the present work under the influence of 3D non-uniform magnetic field.Four distinct Direct Current finite solenoids are used as the source of magnetic field and different Re magnitudes are examined using various flow analysis techniques to understand the complex characteristics.Overall, it is found that, the existence of magnetic field gradients in both radial and axial direction leads to the presence of Kelvin Force Density which significantly alter the flow distribution within the ferrofluid.
Globally, the intricate and distinct gradients of each solenoid configuration are accountable for the occurrence of "toroidal" vortex and "thin cylindrical ring" vortex for distantly spaced and closely spaced solenoid configurations respectively.The critical point theory analysis reveals that, ferrofluid near both the radial extremes demonstrate the "Node/Saddle" dominated topologies, whereas a clear preference is observed for the "Foci" topology at the central radial region (2ρ/D h = 0.5).In addition, existence of two different types of fluid deformations namely pure strain motion and solid body rotation within the same ferrofluid domain, indicates the varying impact of the inertial forces and the magnetic body forces.
Locally, the streamlines demonstrate the occurrence of additional fluid rotations for closely spaced solenoid arrangements along with the existing recirculation zones.The polar distribution identifies the presence of "Cardioid" and "Subcardioid" patterns for the radial KFD component along with the "Omni-directional" patterns for the axial Kelvin Force Density.This profound observation can be utilized in many practical applications where the precise flow monitoring is required along the angular direction.Also, the existence of substantial pressure gradients in the radial direction confirms its equivalent contribution in modifying the flow field.
Finally, based on the flow insights obtained from the present work, it is proposed that a thin finite solenoid with the desired length to diameter ratio can be effectively used to control and manipulate the ferrohydrodynamics based transport devices.

Appendix B. Magnetic field distribution of a solenoid
The OpenFOAM environment does not have an inbuilt way to solve some of the advanced mathematical functions such as complete elliptic integrals, hypergeometric functions, and the Legendre polynomials etc.Thus, its equivalent formulations are solved in terms of the Bessel functions, to compute the magnetic field distributions.These obtained results are then compared with the inaugural work of Callaghan and Maslen (1960) as mentioned in Section 3.1.The distributions are defined for an azimuthally symmetric solenoid of length L s and radius a (D s /2) which consist of N number of turns, and through which a fixed magnitude of current I flows continuously.Both the numerical computations of present study (Fig. 5) are carried out for an iron free, thin finite solenoid using cylindrical coordinates system (ρ, θ, z).
For MATLAB calculations, the field distributions (Eq.B.1-B.5) are provided using complete elliptic integrals based on the work done by Callaghan and Maslen (1960) and Cebron (2020) in which the radial magnetic field variation (B ρ ) is given as Here, K and E refers to the complete elliptic integral of the first kind and second kind respectively, where, Whereas, the axial magnetic field (B z ) is defined based on the maximum value of a non-dimensional parameter (ζ + ), which is described for a combination of radial as well as axial location.
For ζ + < 1, For OpenFOAM computations, both the radial and the axial field variations of a solenoid are solved in terms of the Bessel functions (Labinac et al., 2006;Papadopoulos et al., 2012)

Fig. 4 .
Fig. 4. Grid distribution (a) Along any x − y plane (b) Along any radial line (c) Along any x -z /y -z plane (d) Along any axial line.

Fig. 5 .Fig. 6 .
Fig. 5. Comparison of magnetic field distribution from the present work with Labinac et al. (2006) in (a) Axial direction for 0 ⩽ z < L s (b) Axial direction for z > L s (c) Radial direction.

Fig. 7 .
Fig. 7. Iso-surface of λ 2 = − 0.01 for different L s /D s at Re = 0,Rm = 1 × 10 8 , and ϕ = 0.05.The colour scale represents the non-dimensionalised axial component of KFD, whereas the black solid lines refers to the corresponding Ls start and Ls end .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .Fig. 9 .
Fig. 8. Iso-surface of λ 2 = − 0.01 for different Re at solenoid configuration S1, Rm = 1 × 10 8 , and ϕ = 0.05.The colour scale represents the non-dimensionalised streamwise vorticity, whereas the black solid lines refers to the corresponding Ls start and Ls end .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Sample representation for the type of fluid deformations.

Fig. 14 .
Fig. 14.Iso-magnetic field lines at x mid-plane for different L s /D s (a) S1 (b) S2 (c) S3 (d) S4.The dashed vertical lines represents the corresponding Ls start and Ls end , and the coloured contour represents the H * z .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4)In addition to the previously defined parameters, Π is the complete elliptic integral of the third kind, and as shown in Eqs.B.6-B.8.k; z) J 0 (kρ) J 1 (ka) (B.7)where f(k; z) is defined as a piece-wise function for different values of z,f (k; z) = { e − k(z− Ls) − e − kz , z ⩾ L s 2 − e − k(Ls − z) − e − kz , 0 ⩽ z < L s (B.8)Here, J n for n ⩾ 0 are the Bessel functions of the first kind.
theoretically studied ferrofluids in an oscillating field, and stated that the field dependent viscosity significantly influence the flow vorticity.A similar configuration is S.Dalvi et al.

Table 1
Different solenoid configurations considered in the present study.
Fig. 1.Schematic of the computational domain.

Table 2
List of constant parameters used in the present study.

Table 4
Details of grid independence study.
a Number of nodes for 0.5D h × πD h × 15D h respectively.bThefirstcell height near both the ends of solenoid (Ls start and Ls end ).cThefirst cell height near the cylinder wall (2ρ/D h = 1).S.Dalvi et al.