Mathematical Modelling of an Elongated Magnetic Droplet in a Rotating Magnetic Field

Dynamics of an elongated droplet under the action of a rotating magnetic field is considered by mathematical modelling. The actual shape of a droplet is obtained by solving the initial-boundary value problem of a nonlinear singularly perturbed partial differential equation (PDE). For the discretization in space the finite difference scheme (FDS) is applied. Time evolution of numerical solutions is obtained with MATLAB by solving a large system of ordinary differential equations (ODE).


Introduction
In the second half of the last century, the physics of liquid magnetic materials has become more and more popular. Ferromagnetic particle-based suspensions have achieved the wide usage within the automobile brake and transmission systems. Ferro-fluids, which are the colloids of the ferromagnetic nanoparticles, are being used to produce the seals, to cool the audio speaker coils, as well in biomedicine. Nowadays there is an increasing interest in various soft materials that are created on the basis of the magnetic nanoparticles and have various possibilities of application in mechanics and biomedicine. For these reasons the intensity of the research in the field of soft magnetic materials is rapidly increasing.
The studies of dynamics of the magnetic liquid droplet in a rotating magnetic field was pioneered in [1]. An intriguing reentrant droplet shape transformation was described, in the course of which the oblate shape assumed by the droplet at small strength of the rotating field gets unstable at the critical *  value of magnetic field strength, exceeding which the droplet takes a worm-like shape. The oblate shape is re-assumed at reaching the second critical value of the field strength, however, with a much larger radius than at small field strength. This experimental observation was in remarkable agreement with the theoretical analysis of the stability of equilibrium shapes of droplet in rotating field carried out by virial method [1] (more details of this analysis may be found in [3]). Later this phenomenon was confirmed for other types of magnetic liquid in [9,10]. In addition to the previous findings in [9,10], it was found that the oblate-prolate transition is supercritical for the magnetic permeability µ less than 11 and subcritical for µ > 11. Similar phenomenon is remarked for reentrant transition prolate-oblate at large field strength, which becomes subcritical at µ > 14.
Apart from the oblate-prolate-oblate reentrant transition there is a lot of complex phenomena that take place above the critical values of the field strength. Due to axisymmetric shape instability at field strength above the critical for prolate-oblate transition, the star-fish configurations of the magnetic fluid droplets are formed [1]. Complex dynamics takes place in the range of field strengths where prolate worm-like droplet shape exists [1]. At intermediate field strengths S-like droplet configuration is formed [11], which may break into three daughter droplets [6]. The spirals may be formed on both ends of the droplet [6]. The phase diagram of the experimentally observed shapes of the droplets in a rotating magnetic field is given in [6].
In [2] the mathematical model for the description of the dynamics of wormlike magnetic liquid droplet shapes based on the slender body approach was proposed. The resulting equation for the tangent angle of the center line of the droplet that contains contributions due to viscous, capillary and magnetic forces was regularized introducing the small term due to intrinsic viscosity of the droplet. As a result the steady shapes formed by the propagating shock waves of tangent angle similar to observed in experiments [6] were obtained by numerical calculations.
Here a MATLAB simulation of magnetic droplet dynamics is carried out. The mathematical model from [2] is approximated by using the finite difference scheme (FDS) and method of lines. The central difference FDS is used to approximate the nonlinear second order and the linear fourth order differential spatial operators. The obtained system of ODEs is solved by using MATLAB tool.

Mathematical Model
In [2], the cross section of the droplet of the length 2L is assumed to be circular with constant radius a. The tangential forces along the droplet are neglected. The shape of the droplet is described by the position of its centerline. The PDE for the tangent angle β = β(l, t) of the center line of the magnetic droplet under the action of capillary, magnetic and viscous forces is derived in the following form: where l ∈ (−L, L), t ∈ (0, t f ), here l, t, t f are the arc length of the droplet, the time and the final time moment, β = ωt − θ is the phase lag with angle θ between the local tangent to the centerline of the droplet and the abciss axis, ω is the angular frequency, M = 2π 2 χ 2 H 2 0 a 2 /(µ + 1) is the magnetic torque, δ = 4πη ln(L/a)+c is the hydrodynamic drag coefficient, η i is the intrinsic viscosity of the droplet, γ is the surface tension, χ, µ are magnetic susceptibility and permeability of magnetic liquid respectively, H 0 = (cos(ωt), sin(ωt)) is the rotating field, η is the viscosity of the surrounding fluid, c is a constant of the order of unity.
A shape of the droplet from the known tangent angle is calculated by integration of the following equations The unknown constants of integration are determined from the condition that mass center of the droplet is motionless. In order to present the equation (2.1) in dimensionless form, the characteristic time scale τ = δL 2 /M is introduced. Arc length is scaled by L. As a result the equation reads is the magnetic Bond number given by the ratio of the magnetic and capillary forces, l ∈ (−1, 1), t ∈ (0, t f ).
The regularization term in (2.3) is added from physical considerations using the analogy between bending stress in elastic and viscous filaments [12]. This gives the estimate of parameter ≈ 10 −4 , which works remarkably well as the numerical experiments show. We note that the function F (β) is not monotonic for values of Bm larger than 0.5 (see Fig. 1).
The equation (2.3) is supplemented with boundary conditions corresponding to the absence of normal forces and torques at the ends of the droplet.
By setting = 0 and shifting l to l + 1, we obtain the following problem: In the absence of the last term in equation (2.3) the problem (2.4) is ill posed.

Solution of the Problem
In [2] the numerical solution of (2.3) is obtained by using an implicit scheme with the spatial derivatives approximated with central differences. The nonlinear term due to the function F (β) is resolved by the Newton iterations at each time step. Numerical method for the reduced problem (2.4) shows instabilities.
The stationary solution β s (l) of the equation (2.3) with the boundary condition β(±1, t) = 0 is given by F (β s ) = 0.5ωτ (1 − l 2 ). For the problem (2.4) we obtain that F (β s ) = 0.5ωτ l(2 − l). The maximal value β m is obtained as the solution of the transcendental equation The angle β as a function of the arc length variable l is discontinuous for ωτ > 2F (β 0 ), where β 0 are the roots of equation F (β) = 0 (they define the local maxima of the function F (β)). The values w c = (ωτ ) 0 = 2F (β 0 ), define the critical frequencies.
We modify the nonlinear function F (u) in the following way (see Fig. 14): This process can be continued. Therefore in the segment − 0.5 arccos(0.5/Bm) are the local maxima of the function F (u). The ends of the segment u 2k−1 , u 2k satisfy the following conditions: The critical frequencies w c are defined by the expression w c = 2F (u 2k−1 ), k = 1, 2, . . . . As an example, for Bm = 1.5 we have: The modified function F (β) is monotonic (0 ≤ F (β) ≤ 2+1/Bm) and from (2.4) we can obtain that for fixed time t the solution β(l, t) is quadratically integrable together with their first order generalized partial derivatives with respect to l or β, i.e. it belongs to the Sobolev space The problems (2.3), (2.4) are solved by the MATLAB, using the method of lines and finite difference method for the approximation of the spatial derivatives.
A full analysis of the ill-posed problem (2.4) is beyond the scope of the present work. Here we note that the discretization by finite differences and regularization of function F introduces some numerical bending elasticity which stabilizes the solution for = 0.

Approximation
We consider the uniform grid in the space l j = jh, j = 0, N , N h =L. Using the central finite differences of second order for approximation of partial derivatives of second and fourth order with respect to l, we approximate differential equation (2.3) and boundary conditions β = ∂ 2 β ∂l 2 = 0 by a semi-discrete problem. We write the initial value problem for the system of nonlinear ordinary differential equations (ODEs) of the first order in the following matrix form The discrete approximation of (2.4) is obtained from (4.1) by taking = 0. The equation (4.1) can be rewritten in the forṁ It is well-known that approximation of the second order derivatives by central differences introduces an artificial numerical diffusion where F j = F (β j , t), β j = β(l j , t).
In a similar way a high-order approximation can be constructed for differential problem (2.3) and (2.4). In the case = 0, the following initial value problem for a system of nonlinear ODEs is obtained: The finite difference approximations can be modified using schemes with the exact spectrum (FDSES) [8]. Then in (4.1) the matrix A is replaced with the matrix P DP , where P = P −1 is the symmetrical orthogonal matrix with elements p j,k = 2 N sin πjk N , j, k = 1, N − 1 and the diagonal matrix D contains the N − 1 eigenvalues d k = (kπ/L) 2 of the differential operator (− ∂ 2 ∂l 2 ). Note, that FDSES are also used in [7] to solve hyperbolic heat conduction problems. It is shown in [7] that such schemes can be more efficient than the classical finite difference schemes [4,5] and they can add a numerical diffusion required to regularize ill-posed problems.
For modified vector-function F (U ) (see Fig. 14) the ODEs of (4.1) are implemented in the following form: The shape of a droplet in the plane x, y is found by numerical integration of equation (2.2) x(l, t) = l 0 cos(β(ξ, t)) dξ,ȳ(l, t) = − l 0 sin(β(ξ, t)) dξ, and the center of mass is computed as In the discrete case the trapezoid formula is used to compute integrals: wherex j =x(l j , t),ȳ j =ȳ(l j , t), x j = x(l j , t), y j = y(l j , t), β j = β(l j , t), j = 1, N ,x 0 =ȳ 0 = 0. Therefore, the mass center of the droplet is motionless.              The shapes of a droplet depend on the time and number of jumps of the tangent angle (see Figs. 7, 11, 13). The S-like shape (Fig. 7) can be characterized by two jumps, but the 8 type shape (Figs. 11, 13) is characterized by four jumps of the tangent angle. The obtained shapes are remarkably similar to shapes observed in real experiments [6].
In Tables 1, 2, 3 we can see the maximal value β m of β and the number of time steps K at t = t f in dependence on the parameters and ωτ . All results are obtained with the Matlab solver "ode15s" (RelTol = 10 −6 ) correspondingly from (4.1) (FDS-O(h 2 )), (4.2) (FDS-O(h 4 )), (4.1) (FDSES). The maximal values of β m or β(1, t f ) and the number of time steps K are decreasing in dependence on the parameter .

Conclusions
It is found that the ill-posed problem for parabolic type PDEs may be efficiently solved by applying the regularization of the diffusion functional and using the MATLAB solver "ode15" to the nonlinear system of ODEs obtained after approximation of the spatial derivatives with finite differences. The transition of simple droplets to more complex droplet shapes is investigated by simulating the propagation of jumps of the tangent angle.
The developed numerical method allows us to resolve the shapes of the ferrofluid droplet in a rotating magnetic field observed experimentally, such accuracy of simulations was not possible by using the conventional method of regularization based on physical considerations.
A simple model of the elongated magnetic liquid droplet rotating synchronously with an applied magnetic field has been proposed. Temporal dy-namics of the droplet shape evolution for the given parameters and the stationary state solutions corresponding to S-like shape are obtained. With the increase of the frequency of a rotating field several shock waves develop. The ends of the droplet become spiralized at high frequencies. At the critical value a jump of the tangent angle appears near the droplet's center and propagates until a new steady shape of the droplet is established. The mathematical model allows us to identify the critical frequencies at which the droplet's shape transitions take place. Shapes obtained numerically are remarkably similar to the ones found in experiments [6].