Numerical Simulation of Magnetic Droplet Dynamics in a Rotating Field

Dynamics and hysteresis of an elongated droplet under the action of a rotating magnetic field is considered for mathematical modelling. The shape of droplet is found by regularization of the ill-posed initial–boundary value problem for nonlinear partial differential equation (PDE). It is shown that two methods of the regularization – introduction of small viscous bending torques and construction of monotonous continuous functions are equivalent. Their connection with the regularization of the ill-posed reverse problems for the parabolic equation of heat conduction is remarked. Spatial discretization is carried out by the finite difference scheme (FDS). Time evolution of numerical solutions is obtained using method of lines for solving a large system of ordinary differential equations (ODE).


Introduction
Dynamics of magnetic droplets in the rotating magnetic field is an exciting field of research. In [1] it was found that magnetic droplets with high magnetic permeability have sequence of shape bifurcations with the increase of the magnetic field including transition from oblate shape to prolate one and back to oblate at higher values of the magnetic field strength. These experimental results are in good agreement with the theoretical calculations by the virial method [1] (for more details see [4]). Later in [9,10] it was shown that if the magnetic permeability of the droplet is enough high these bifurcations are subcritical and therefore hysteresis phenomena should occur. Behavior of the droplet at intermediate values of the rotating field strength were investigated in [7] where further bifurcations of elongated droplets were observed, -formation of shapes with large curvature in definite places of the droplet corresponding to the discontinuities of the tangent of the center line of elongated droplet. This behavior was theoretically reproduced by the model proposed in [2], which consider the dynamics of the droplet due to the surface tension forces and torques due to the rotating magnetic field. The ill-posed nonlinear PDE for the tangent angle was derived in [2]. Its numerical solution was regularized by introducing the small viscous torques due to the bending of the droplet [13]. In our paper [3] it is shown that the regularization of the ill-posed problem may be achieved by introducing modified continuous monotonous function with discontinuous first derivative which allow us to consider the path with increasing frequencies of the rotating field. Here by applying this technique two modified functions are constructed which allow us to consider the paths with increasing and decreasing frequencies of the rotating field. As a result the existence of the multiple stationary states of the droplet in definite ranges of the frequency of rotating field and therefore hysteresis phenomena are predicted.
In the paper [5], the regularization techniques for the backward in time nonlinear parabolic problem are considered.

Mathematical Model
In [2] the cross section of the droplet with 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 nonlinear PDE of parabolic type 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 [2,3]: where the parameters l, η i , t, t f are the arc length, the intrinsic viscosity of the droplet and time with the final time moment t f , β = ωt − θ is the phase lag with angle θ between the local tangent to the centerline of the droplet and the abscissa axis (ω is the angular frequency of a rotating field), M = is the magnetic torque, δ = 4πη ln(L/a)+c is the hydrodynamic drag coefficient, η i , γ are the intrinsic viscosity and the surface tension of the droplet, χ, µ 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.
Shape of the droplet depending on the time in the frame of rotating field is calculated from the known tangent angle β by integration of the set of ordinary differential equations (ODEs) 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. We have following PDE: where F (β) = 1 Bm β + sin(2β) is the nonlinear function, = 3πa 4 ηi 4δL 4 is a small coefficient (of order 10 −4 ) for the regularization, Bm = is the magnetic Bond number given by the ratio of the magnetic and capillary forces, The regularization term in (2.3) with parameter is added from physical considerations using the analogy between bending stress in elastic and viscous filaments [13]. The function F (β) is not monotonic for Bm values larger than 0.5 (see Fig. 1 for Bm = 1.5). The equation (2.3) is supplemented by boundary conditions corresponding to the absence of normal forces and torques at the ends of the droplet: β(0, t) = β(2, t) = ∂ 2 β(0,t) The initial condition is defined by β(l, 0) = Θ 0 (l), l ∈ [0, 2], where Θ 0 (l) is continuously differentiable function with Θ 0 (0) = Θ 0 (2) = 0. The last term in the equations (2.1), (2.3) is used for the regularization of the numerical approximations.
By setting = 0 we obtain the following ill-posed problem: In the numerical experiments by changing the frequencies ωτ we can take Θ 0 (l) equal to the stationary solution obtained at the previous assigned value of the frequency.
The second method of regularization of the ill-posed problem (2.4) consists in construction of monotonous continuous functions. For the path with increasing frequency ωτ of rotating field the modified (direct) function F (β) = F (u) is constructed as follows [3]: , where u 1 = π 2 − 0.5 arccos(0.5/Bm), is the first local maxima of function F (u) or the solution of the equations More generally, we can calculate the maxima of the function and the coordinates of the intersection points. Therefore in the segment [u 2k−1 , u 2k ], k = 1, 2, . . . the function F (u) is replaced with line segment F (u) = F (u 2k−1 ) = F 2k−1 , where u 2k−1 = (2k−1)π 2 −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 following conditions: . This function is shown in Fig. 2 for Bm = 1.5.
For the path with decreasing frequency of the rotating field the modified (reverse) function F (β) = f (u) is defined as follows: is the value of the local minimum of function f (u) or the solution of the equations f (u) = 0 in the segment [u 3 , u 4 ], Fig. 3.
The stationary shapes with 3 jumps constructed according to the modified functions are shown in Fig. 4 (direct function) and Fig. 5 (reverse function) for Bm = 1.5, ωτ = 15. It is interesting to remark that the curvature of the shapes has opposite signs for the cases of direct and reverse functions. From For the stationary solutions depending on the frequency we can obtain one or two solutions. In the segments ωτ ∈ [2F 2k−1 , 2f 2k−1 ], k = 2, 3, . . . we have two stationary solutions, but in the segments ωτ ∈ (2F 2k−1 , 2f 2k+1 ), k = 2, 3, . . . the solution is unique (F (u) ≥ 0). An example, at Bm The problem (2.3), which is regularized by , and problem (2.4), which is regularized by the modified functions, are solved by the MATLAB "solver ode15s" with relative error 10 −7 (RelTol = 10 −7 ), using the method of lines and finite difference method for the approximation of spatial derivatives. We consider the uniform grid in the space l j = jh, j = 0, N , N h = 2. The finite differences of second order approximation for partial derivatives of second and fourth order with respect to l are applied. As a result the initial value problem for the system of nonlinear ODEs of the first order is obtained: For the difference scheme with exact spectrum (FDSES) [3] the matrix A is replaced with the matrix W DW , where W = W −1 is the symmetrical orthogonal matrix with elements w j,k = 2 N sin πjk N , j, k = 1, N − 1 and the diagonal matrix D contains the N − 1 and eigenvalues λ k = ( kπ L ) 2 of the differential operator (− ∂ 2 ∂l 2 ). For FDS the diagonal matrix D has eigenvalues µ k = 4 h 2 sin 2 ( πk 2N ). Using the matrix form A + h 2 12 B for the fourth order approximation of the diffusion operator [12], we can obtain the following problem: It should be remarked that the regularization of ill-posed backward in time linear homogeneous heat transfer equation (2.3) by term ∂ 4 β ∂ 4 l or BU (t) for ODEs system (3.1) is considered in the book by Lattes and Lions [8].
Inverse problems for partial differential equations and their methods of regularization are considered in the book by V.Isakov [6].
By using the Lattes and Lions regularization we will consider the following two initial value problem: and Then we get the problem: from which the following integral identity follows trivially: where g(u) = 2 3 + 2 cos(2u), − 4 3 ≤ g(u) ≤ 8 3 , g (u) = −4 sin(2u). In order to investigate the solvability of the corresponding initial-value problem in the Sobolev space 0 W 2 2 , to prove the existence of the weak solution and to obtain a priori estimations at fixed time t, we need to determine the parameter from the following inequality: Since φ is arbitrary, then for fixed t the following nonlinear differential equation is obtained This equation is solved numerically by using Matlab solver "bvp4c", the following 5 boundary conditions are used The last condition is used to find κ for the fixed value of parameter k 0 . The maximal value of κ = 0.0143 is obtained for k 0 = 1.3. Then from (3.6) and using Hölder's inequality Using the Friedrich inequality [11] u 2 ≤ 4 π 2 u x 2 we obtain the inequality from which it follows that t .
In [8] the equation (3.8) is replaced by the regularized high-order equation Then a k (t) = exp(g(λ k − λ 2 k )(T −t))a T k and the series for u(l, t) is convergent in L 2 (0, 2) for any u T and any t < T . Moreover, when goes to 0, the regularized solutions are converging to the exact solution u of the initial problem (3.8) [6].
The second method for regularization of (2.4) is given by the initial value problem (3.1) with = 0 but with the modified functions F , f. The shape of a droplet in the plane x, y is found by numerical integration of the equation (l, t) dl.
. Therefore, the mass center of the droplet is motionless. In the discrete case the trapezoid formula is used.

Numerical Results
Constructed direct and reverse functions allow us to calculate the dynamics of shapes corresponding to the path with increasing the frequency of rotating field starting from straight configuration and the reverse path starting from the deformed shape calculated by using direct function. The main numerical simulations are carried out by N = 100 and integrating the system of ODEs

Numerical results for different frequencies using the direct path
The dynamics of the tangent angle and the corresponding shapes formed in direct path obtained for Θ 0 (l) = 0 for different frequencies of rotating field by Bm = 1.5, ωτ = 5, 8, 12, 24 are shown in Figs. 6-13. Here the function β(l, t) with 1, 2, 3 and 5 jumps at different time-moments t is shown (the stationary solution is obtained by t f = 6). We note the formation of highly spiral -like shapes at large frequencies of the rotating field (see Figs. 11,13).
In [3], by using the Matlab solver "ode15s" we have solved the same problem by taking = 0. Now we compare more models, including different regularization techniques. The numerical results for parameters ωτ = 8, N = 100 and different systems of ODEs (3.1), (3.2), (3.3), (3.4) with = 10 −4 and = 0 are presented in Table 1. Here we analyze the dependence of the maximal value β m of β and the number of time steps K at t = t f = 6 on the different integration and regularization methods and values of parameter ε. The results are obtained with the Matlab solver "ode15s" (RelTol = 10 −7 ). For > 0, we have    observed that the computing process in time is smooth (see also Figs. 16, 17), but for = 0 we have observed the oscillations in time (these cases are denoted by * in Table 1, also see Figs. 14, 15).

Numerical results for different frequencies using the direct and reverse paths
The main results of this paper are shown in Figs. 26-31. In this case the shape obtained by reverse function starting from the stationary configuration with t f = 6, which is obtained by using the direct function at higher frequency, is different from the configuration obtained at the same frequency by the direct path. So shapes obtained at ωτ = 2 by direct and reverse paths are different (see Fig. 26). Here we plot two stationary solutions at frequency ωτ = 2 obtained at t f = 6 in the following way: 1) We use the direct path (function F ) with Θ 0 (l) = 0 and obtain the stationary solutions by ωτ = 2 and ωτ = 5; 2) We use the reverse path (function f ) with Θ 0 (l) equal to the stationary solution obtained with the direct function for ωτ = 5 and obtain the stationary solution for ωτ = 2.
We will denote this scenario as: from 0 to 2 and to 5, then from 5 to 2 or ωτ = 0 ⇒ 2 ⇒ 5 ∪ 5 ⇒ 2. The same is valid for ωτ = 6, ωτ = 0 ⇒ 6 ⇒ 8 ∪ 8 ⇒ 6 ( Fig. 28) and ωτ = 10, ωτ = 0 ⇒ 10 ⇒ 12 ∪ 12 ⇒ 10 (Fig. 30). The different tangent angle  Thus the considered model predicts multiple stationary states of the droplet in definite ranges of the frequency of rotating field. It would be interesting to confirm this prediction in experiments. Here we should remark that available experiments [7] indeed show the formation of the shapes with discontinuity of tangent angle which breaks at places with large curvature. The breaking phenomenon is not described by the present model.

Conclusions
The ill-posed problem for non-linear parabolic PDE may be regularized by introducing viscous bending torques and construction of monotonous continuous functions. By numerical simulation it is found that both approaches are equivalent. It is shown that increasing the accuracy of approximation of the discrete problem allows to decrease the number of time steps to obtain the same results.
Two monotonous functions corresponding to the paths with increasing and decreasing frequencies may be constructed and two different droplet shapes exist in definite ranges of the frequency of the rotating field. This leads to hysteresis phenomena in the droplet shape transformations at change of the frequency of rotating field. This prediction is a challenge for experimental verification.