ADI Approach to the Particle Diﬀusion Problem in Magnetic Fluids ∗

. The present study is devoted to the development of an ADI approach to simulate two-dimensional time-dependent diﬀusion process of ferromagnetic particles in magnetic ﬂuids. Speciﬁc features of the problem are the Neumann boundary conditions. We construct an ADI scheme of formally second order accuracy approximation in time and space. It is proved that the scheme is absolutely stable and it has the accuracy in the energy norm of the second order in time and the order 32 in space. The numerical results of a test problem indicate that the convergence rate in space is of the second order as well.


Introduction
In recent years, papers related to a study of the diffusion process of ferromagnetic particles take a considerable place in the hydromechanics of magnetic fluids.As known, a magnetic fluid is a stable colloidal suspension of Brownian suspended particles of ferromagnetic in a non-magnetic carrier liquid.Under the action of a nonuniform magnetic field, the process of particle diffusion with respect to the carrier liquid takes place, resulting in a redistribution of the fluid magnetization and exerting thereby a significant influence on the hydrostatics of magnetic fluids.In [8], the exact solution of the steady-state particle concentration problem has been given for any closed computational domain in any space dimension.Using this solution, two ferrohydrostatic problems on the influence of the particle diffusion have been solved numerically: on axisymmetric equilibrium shapes of a free magnetic-fluid surface in the magnetic field of a cylindrical conductor with direct current [3] and on the instability of a horizontal magnetic-fluid layer in a uniform magnetic field [7].In [1,2] twodimensional time-dependent problems on steadying the particle concentration distribution with time in a rectangular magnetic-fluid domain are solved by means of a finite-difference scheme of alternating direction (ADI) type.
The present study is devoted to the development of an ADI approach to time-dependent particle diffusion problems, specific features of which are the boundary conditions of the Neumann type.It should be noted that the problem of constructing appropriate boundary conditions for splitting difference schemes is not easy even in the case of Dirichlet conditions, and requires a very special analysis [4,6].Before, efforts have been undertaken to build ADI schemes for parabolic-type PDE with Neumann conditions, e.g. in [5] for the heat equation, in [14] for the convection-diffusion equations, in [13] for the phase field model equations.Increase of the approximation error order at boundary nodes has been achieved by means of the usage of the given differential equation.We consider a Neumann problem for the two-dimensional particle concentration equation modified by a special change of variables.The ADI technique is applied to solve the problem for the first time.Notice that in the case of Dirichlet boundary conditions, the ADI scheme for the modified concentration equation is studied in [12].

Mathematical Model
The magnetic particle mass transfer in a magnetic fluid (the mass conservation law for Brownian magnetic particles) can be described in dimensionless form (see [1,2,8]) by where C is the relative volume concentration of the particles, t the dimensionless time, ξ = ξ(x) the dimensionless intensity of magnetic field, x the space coordinates, µ 0 = 4π × 10 −7 H/m the magnetic constant, m the magnetic moment of a particle, H c a characteristic magnetic-field intensity in the fluid volume (e.g. a maximal intensity value), k = 1.3806568 × 10 −23 J/K the Boltzmann constant, T the temperature of a particle, L(α) = coth α − 1/α the Langevin function, A a dimensionless magnetic parameter.Equation (2.1) is supplemented by the condition describing the impermeability of boundaries by particles, i.e. the particle flux through boundaries is equal to zero: which is a boundary condition of the Robin type.
As the initial condition, we consider a uniform concentration of particles where V denotes the fluid domain, |V | its volume, and the coefficient b is a positive function of spatial variables.
From the numerical modelling point of view, it is not acceptable that the normal derivative of magnetic intensity in condition (2.2) is a sign-changing function, because it can lead to computational instability.In order to improve the mathematical model, we introduce a change of variables In this paper we are interested in the 2D particle concentration problem posed in a rectangular domain Then, the problem with respect to the function u(x, t) is given in the following form where the Neumann boundary condition is given at the boundary Γ of the rectangle Ω.

ADI Scheme
We consider a more general problem where ) and (3.1) coincide.
On the rectangle Ω we introduce the uniform mesh with the step sizes h 1 and h 2 in x 1 and x 2 directions, respectively.
For the approximation of derivatives by finite differences and for mesh function values, we shall use the notation introduced in [11].For example, u xα and u xα define the forward and backward difference quotient of u with respect to the variable x α , b (±0.51) = b(x 1 ± 0.5h 1 , x 2 ), b (±0.52) = b(x 1 , x 2 ± 0.5h 2 ), ϕ j+1/2 = ϕ(x, t + 0.5τ ) where τ is the mesh step by time.
Thus, problem (3.1) is approximated on the uniform space-time mesh by using standard symmetrical finite differences in x 1 and x 2 directions.The boundary conditions are approximated with a local truncation error O(τ 2 + |h| 2 ), |h| = max(h 1 , h 2 ), assuming that the differential equation (3.1) is valid at the boundary nodes.For example, the boundary condition at x 1 = 0 is approximated as .
Then, we represent the approximated boundary condition as follows extracting the time derivative in explicit form, similarly to the representation of the differential equation at inner nodes.The other boundaries and angular points are handled in the same way.For example, at the lower left angular point we get for i 1 = 0, i 2 = 0: The main objective of this paper is the construction of an ADI scheme for problem (3.1).We consider the following ADI scheme b(x) where for α = 1, 2, (3.4) Math.Model.Anal., 16(1):62-71, 2011.
It is easy to check that scheme (3.2)-(3.4)has a local truncation error of order O τ 2 + |h| 2 on the mesh ω including all boundary nodes.Since it is of the form (3.2) not only at internal mesh nodes but also at boundary nodes, the three-point Thomas algorithm is applied not only along internal lines of the mesh as in the case of ADI scheme for Dirichlet boundary conditions [12], but also at the boundaries of the mesh domain.Note that the Thomas algorithm for scheme (3.2)-(3.4) is absolutely stable.

Stability of the ADI Scheme
In order to prove stability of the scheme (3.2)-(3.4),we use the following lemma given in [12]: where It is not difficult to verify directly that A α = A * α ≥ 0, B = B * > 0, i.e. the difference operators A 1 and A 2 are self-adjoint and nonnegative definite, and the operator B is self-adjoint and positive definite in the space of mesh functions H defined on the mesh ω.For example, let us prove the non-negativity of the operator A 1 .We define an inner product in H as a discrete analog of the inner product in the Hilbert space L 2 (Ω): In order to prove that [A 1 y, y] ≥ 0 for all y ∈ H, it is sufficient to show that [A 1 y, y] i2 ≥ 0 for any 0 ≤ i 2 ≤ N 2 .Indeed, using the first Green difference formula [11] we obtain As a norm we choose y (1) := (B + 0.5τ A 2 )y D which is used in [11,12].Applying estimates (4.1) and (4.2), we get the chain of inequalities
We investigate the convergence of the scheme (3.2)-(3.4)for ϕ = 0 which corresponds to the diffusion problem (2.5).Let z = y − u be the error of the difference solution.It is obvious that for ϕ = 0 the mesh function z satisfies the problem where The problem (5.2) is of the form (5.1).In view of inequality (4.3) and using the assumption z 0 = 0, we get the following estimate At internal nodes, the mesh function ψ represents the approximation error of the schemes (3.2) and (5.1) for ϕ = 0, i.e. ψ = O τ 2 + |h| 2 .At the left boundary, it holds Similar relations can be shown on the other parts of the boundary.At the vertices of the domain, we have ψ = O τ 2 + |h| .Then, we conclude From (5.3) it follows that

Test
The convergence order obtained theoretically does not coincide with the approximation order of the constructed ADI scheme, thus it is reasonable to check this fact by numerical experiments.As a test problem, the particle diffusion problem in a square cavity under the action of a magnetic field of a permanent magnet with a hyperbola-shaped polar head is chosen (see Fig. 1).The diffusion process is simulated using the solution u(x, t) of the problem (2.5), where C(x, t) = b(x)u(x, t), b = sinh(Aξ) (Aξ).The magnetic field intensity is described by the formula [9,10] .
In the steady-state, the error z(x) of the numerical solution C(x) is computed as z = C − C exact , where C exact is the exact solution of the steady-state problem which, according to formula (2.4), can be written as In order to determine the order of convergence of the difference solution to the exact solution as N → ∞, the dependence has been considered for different meshes, and the unknowns a and p are determined for each pair of neighbouring values of N .
The results of the test are presented in Tables 1 and 2. They show that for the given test problem scheme (3.2)-(3.4)demonstrates the second order accuracy in space.
We also have tested the constructed ADI scheme on two model problems of generalized form (3.1) with known exact solutions, which are not solutions of the concentration problems.Our computations confirmed that the scheme demonstrates the second order accuracy in space and time.

Conclusions
In this paper, an ADI scheme for the 2-D time-dependent magnetic particle diffusion problem defined in a rectangular magnetic-fluid domain with Neumann boundary conditions, is constructed, analyzed and tested.The theoretical analysis shows that the scheme is absolutely stable and converges with the second order in time and 3  2 order in space.For test problems, we found numerically that the scheme demonstrates the second order in space as well.We hope to prove this fact analytically in near future.

Figure 1 .
Figure 1.Illustration of the problem statement.

Table 1 .
Convergence rate in the energy norm.

Table 2 .
Convergence rate in the maximum norm.