Numerical Study of Diffusion of Interacting Particles in a Magnetic Fluid Layer

Magnetic fluids are stable colloidal suspensions of ferromagnetic nano-particles (of size 10-20 nm) in a nonmagnetic liquid carrier. An initially uniform particle distribution in the carrier becomes spatially inhomogeneous in nonuniform magnetic fields. Motion of particles in magnetic fluids under the action of magnetic fields is of particular interest for contemporary mathematical and numerical modelling in ferrohydrodynamics.


Introduction
Magnetic fluids are stable colloidal suspensions of ferromagnetic nano-particles (of size 10-20 nm) in a nonmagnetic liquid carrier. An initially uniform particle distribution in the carrier becomes spatially inhomogeneous in nonuniform magnetic fields. Motion of particles in magnetic fluids under the action of magnetic fields is of particular interest for contemporary mathematical and numerical modelling in ferrohydrodynamics.
The most theoretical models for the diffusion process in magnetic fluids assume no interaction between particles (Bashtovoi et al., 2007;Lavrova et al., 2010;Polevikov & Tobiska, 2008;, which is valid for dilute fluids only. This assumption allows to construct an explicit dependence between equilibrium particle concentration and the magnetic field distribution, simplifying significantly the modelling. In case of concentrated magnetic fluids another theoretical model should be considered. Recently, a dynamic mass transfer equation for describing diffusion of interacting ferromagnetic particles in magnetic fluids was derived (Pshenichnikov et al., 2011). In this paper it is mentioned that "...In the case of high particle concentrations, the magnetic and diffusion problems are strictly interrelated, and the concentration profile depends markedly on steric, magnetodipole, and hydrodynamic interparticle interactions, whose counting is a problem of great concern..." The present study is devoted to the classical problem of ferrohydrostatics on stability (known as the normal field instability or the Rosensweig instability) of a horizontal semi-infinite layer of a magnetic fluid under the influence of gravity and a uniform magnetic field normal to the plane free surface of the layer (Rosensweig, 1998). A periodic peak-shaped structure is formed on the fluid surface when the applied magnetic field exceeds a critical value. This phenomenon was observed first experimentally (Cowley & Rosensweig, 1967).
A number of papers are devoted to numerical investigations of the Rosensweig instability. The numerical results concern equilibrium states of the ferrofluid layer in (Aristidopoulou et al., 1996;Bashtovoi et al., 2002;Boudouvis et al., 1987;Gollwitzer et al., 2007;2009;Lavrova et al., 2003;2010) and analyze dynamical properties of the ferrofluid 9 www.intechopen.com behavior in (Knieling et al., 2007;Matthies & Tobiska, 2005). A non-uniform equilibrium distribution of particles in the ferrofluid layer is computed in (Lavrova et al., 2010) for the first time.
In order to reach the equilibrium between concentration and the magnetic field, quite a long time is needed. The concentration remains almost constant for much shorter time scales. That is why, the validity of the results, mentioned in the previous paragraph, will not be abolished by the preset contribution. The aging process of the Rosensweig instability was experimentally studied over the long time (up to 50 days) in (Sudo et al., 2006). The effect of evaporation and non-evaporation of magnetic fluid on the pattern aging were examined. It was found that the interfacial spike pattern of magnetic fluids changes with time dramatically. Namely, the cell pattern gradually bifurcates at the constant magnetic field, and the number of spikes increases with time.
The particles are in Brownian motion inside the ferrofluid layer, when no magnetic field is applied, and the particle concentration is constant over the fluid volume. This is correct under an assumption that the gravity force has a negligible influence to the diffusion of particle. When the applied field is switched on but the field intensity is too weak to perturb the plane surface, then the magnetic field inside the layer remains constant and the particle concentration is constant, as a consequence. The situation changes when the applied field intensity is strong enough to perturb the free surface. A nonuniformity of the magnetic field inside the fluid causes a redistribution of the particles. This is due to interactions between field and particles and interparticle interactions. An interaction between particles and the magnetic field of the fluid was taken into account for the modelling of the Rosensweig instability in our previous research in (Lavrova et al., 2010). This interaction plays a dominant role in dilute magnetic fluids. The main objective of the contribution is the extension to the case of interacting particles, which should be taken into account for concentrated magnetic fluids.
Mathematical model of the coupled problem consists of the magnetostatic subproblem, the concentration subproblem and the free-surface subproblem. The concentration subproblem is based on a recently developed mass transfer equation for describing diffusion of interacting ferromagnetic particles in magnetic fluids (Pshenichnikov et al., 2011). Three subproblems are strictly interrelated to each other and should be solved simultaneously to resolve equilibrium states of the system. An iterative decoupling strategy is applied for solving the coupled system of equations. The finite-element method is used for discretization of the magnetostatic subproblem in a fixed domain. The Newton method is applied to find an element-wise distribution of the concentration at the finite-element mesh. The finite-difference approach is used for the free-surface subproblem. Numerical results of three models (model 1nonuniform particle distribution without particle interaction, model 2 -nonuniform particle distribution with particle interaction, model 3 -uniform particle distribution) are compared.
The effect of particle interaction shows a considerable influence on behavior of the ferrofluid layer in a uniform applied magnetic field.

Mathematical model
We consider a semi-infinite ferrofluid layer with a horizontal plane free surface bounded from above by a nonmagnetic gas (air). The unperturbed free surface is defined by equation z = 0, whereas the fluid corresponds to a region z < 0. The system is regarded under the action of gravity g =( 0, 0, −g) and a uniform magnetic field H 0 =( 0, 0, H 0 ) normal to the plane free surface of the layer. We consider a single peak in the surface pattern with a cell Ω cell and a free surface Γ. The problem will be formulated in a cylinder Ω cell × (−∞, +∞), see Fig. 1 for the case of a hexagonal cell. The mathematical model for a non-uniform equilibrium distribution of ferromagnetic particles in a magnetic fluid with a free surface leads to a coupled problem formulation consisting of three subproblems. The first subproblem describes the magnetic field structure inside the fluid and in the surrounding air by the Maxwell's equations. The second subproblem concerns the diffusion of particles in the bulk of the fluid as a steady-state concentration problem. Finally, the third subproblem is given by the generalized Young-Laplace equation for equilibrium free-surface shapes of the fluid-air interface.
The Maxwell's equations inside the magnetic fluid and in the air are ∇×H = 0, ∇·B = 0i n Ω cell × (−∞, +∞),( 1 ) see e.g. (Rosensweig, 1998). Here H denotes the magnetic field, B = µ 0 (H + M) the magnetic induction, µ 0 = 4π × 10 −7 H/m the permeability of vacuum. The magnetization vector M is parallel to the field vector and follows a magnetization law M = M(H, C) dependent on the field intensity H and a particle concentration C. An equilibrium magnetization of magnetic fluids with account for interparticle interaction is derived in (Ivanov & Kuznetsova, 2001) in the framework of the modified model of the effective field

185
Numerical Study of Diffusion of Interacting Particles in a Magnetic Fluid Layer

www.intechopen.com
Here M s is the saturation magnetization, C the volumetric concentration and C 0 = 1 Cdx for the fluid domain Ω f with the volume |Ω f |. L(ξ)=coth (ξ) − 1/ξ is the Langevin function, γ the Langevin parameter, χ L the initial susceptibility of the Langevin magnetization, H e the effective field. The magnetization of air equals to zero. The magnetic field satisfies continuity conditions at the interface Γ between ferrofluid and air, see (Rosensweig, 1998) Here [·] denotes a jump over the interface, τ k and n are tangential and normal vectors to the interface. A symmetry condition is specified at the side of a cylinder domain The uniform applied field H 0 is perturbed only in a neighborhood of the interface Γ.T h a ti s why we introduce asymptotic boundaries z = ±δ, far enough from the interface, and specify there a uniform magnetic field where δ > 0. The intensity of the applied field H 0 presents a control parameter of the model, whereas the field H 1 0 =( 0, 0, H 1 0 ) will be computed from the second transmission condition (3), satisfied at the flat interface z = 0. The Maxwell's equations (1) together with conditions (3)-(5) present the first subproblem of the mathematical model.
The second subproblem describes a magnetophoresis process, i.e. the diffusion of ferromagnetic particles in the magnetic fluid under the action of a nonuniform magnetic field. A dynamic mass transfer equation for describing diffusion of interacting ferromagnetic particles in magnetic fluids was derived in (Pshenichnikov et al., 2011) in Ω f , t > 0. The equation is presented with an assumption that the gravity force has a negligible influence to the diffusion of particles. The constant D 0 denotes Einstein's value of the diffusion coefficient for dilute solutions, is the relative mobility of particles in the magnetic fluid. A function specifies the contribution of a magnetodipole interaction to the free energy of the dipolar hard sphere (Pshenichnikov et al., 2011). Here λ is the dipolar coupling constant or the aggregation parameter, estimating the intensity of the magnetodipole interaction in comparison with thermal energy. The modified effective field model, which is used to describe the equilibrium magnetization (2), is applicable for λ ≤ 2 (Pshenichnikov & Lebedev, 2004). We take λ = 1in our model and get The static distribution of particles in the cavity is obtained by equating the full particle flux to zero. According to (Pshenichnikov et al., 2011), it gives To fix the constant c c the condition of conservation for the concentration over the fluid domain will be used. Let us substitute the function G(1, C) to equation (6), then where R(C) is a rational function of the concentration The second subproblem of the mathematical model is presented by equation (8) and condition (7) to fix the constant c c .

Remark 2.1. Equation (8) can be reformulated as
Such a representation is similar in form to the solution for the concentration in dilute approximation (Polevikov & Tobiska, 2008)

Thus, for concentrated fluids the explicit dependence of C on H (9) is replaced by the implicit dependence (8).
The third subproblem of the model defines a shape of the interface Γ between the magnetic fluid and the air. Equilibrium shapes of a free magnetic-fluid surface are described by the generalized Young-Laplace equation, which presents the force balance at the fluid-air interface Here σ is the surface tension coefficient, K the sum of principal curvatures, p 0 the pressure in the air and p is the fluid pressure. The equation of hydrostatics for magnetic fluids is MdH , see e.g. in (Rosensweig, 1998), where ρ is the fluid density and e z =( 0, 0, 1). This equation allows us to express in explicit form the fluid pressure as where c f is an integration constant. The equation (10) for the fluid pressure (11) is where The constant c f will be determined by integrating equation (12) over the surface Γ,s e e Section 3.3 for details.
A solution of the magnetostatic subproblem (1)-(5) depends on the particle concentration and on the shape of the fluid-air interface. The concentration distribution is a solution of the nonlinear equation (8) and depends on the magnetic field configuration inside the ferrofluid. The fluid-air interface satisfies equation (12), which depends on the magnetic field and the concentration at the free surface of the ferrofluid. Three subproblems are strictly interrelated to each other and should be solved simultaneously to resolve equilibrium states of the system.

Numerical methods and tools
Experiments show that the surface pattern of the Rosensweig instability is presented by a hexagonal or square array of spikes, see e.g. (Gollwitzer et al., 2006). For sake of simplicity, we assume that the cell has a circular shape of radius a inscribed to a hexagonal cell of the pattern, see Fig. 2. It allows us to study axisymmetric solutions of the presented model in a two-dimensional geometry of cylindrical coordinates. Axisymmetric solutions on a circular cell were compared with three-dimensional ones on a hexagonal cell in (Lavrova et al., 2003). It was found that the profile shapes of both models nearly completely coincides. The different geometry of the cell results in a 10 % smaller amplitude of an axisymmetric peak in comparison with a three-dimensional solution.

Magnetostatic subproblem
Let us consider the magnetostatic subproblem (1)-(5) in a domain with the fixed interface and assume that the spacial distribution of the particle concentration is given. We express the magnetic field in terms of a scalar magnetic potential φ as H = ∇φ and reformulate the magnetostatic subproblem (1)-(5). The first Maxwell's equation (1) is exactly satisfied, whereas the second one takes form of an elliptic partial differential equation with jumping coefficients, Here Ω ax := Ω 1 ∪ Ω 2 is a meridional cross-section of the 3D-domain Ω,whereΩ 1 corresponds to the fluid and Ω 2 to the air. Differential operators are though in cylindrical coordinates with an assumption of an axial symmetry for the potential. The magnetostatic problem is closed by a set of conditions The boundary Γ D consists of the top and bottom boundaries of the rectangular domain Ω ax = (0, a) × (−δ, δ) and Γ N = ∂Ω ax \ Γ D . We have from condition (5) that Remark 3.1. The constant H 1 0 is defined from the second transmission condition (14), satisfied for the undisturbed interface z = 0. The potential is given in this case as We get The nonlinear equation (17) is solved by the Newton method. Starting from the value H 0 ,themethod converges for 3-5 iterations for the relative error 10 −10 .
Due to the problem reformulation in cylindrical coordinates and the assumption of axial symmetry, a corresponding variational problem is formulated in weighted Sobolev spaces: A structured triangular mesh is used for discretization. A schematic representation of the mesh is shown in Fig. 3. An algorithm of bilinear interpolation is applied for the mesh construction. This algorithm defines every interior grid point (filled circular marker in Fig.  3) by interpolation of eight boundary points (empty circular markers) where parameters s, t ∈ (0, 1) and indices represent the quarter directions of boundary points (N -north, S -south, W -west, E -east), relative to the inner point x.T o p r o d u c e a n interface-fitted mesh, the algorithm is applied in Ω 1 and Ω 2 separately, based on the pointwise representation of the interface, as a solution of the free-surface subproblem, and a uniform distribution of grid points at the boundary sides of Ω 1 and Ω 2 , see Fig. 3. A finite element mesh is reconstructed every time, when the interface is changed. By construction all meshes have the same topology. It allows to define an initial approximation of the potential at the new mesh as a solution of the magnetostatic problem at the old mesh without any interpolation.
The variational problem (18) is discretized by continuous piecewise linear functions on triangles for the given concentration C. Nonlinearities in the discrete equations are treated by a fixed-point iteration One Gauss-Seidel step is applied to the resulting system of linear equations in each iteration. The iterative process (19) continues until the relative a-posteriori error will be smaller than ǫ p (generally 10 −5 ) Here N p denotes the number of unknowns and φ i = φ h (ξ i ) is the nodal value of the potential at the grid point ξ i . The iterative process needs 5-10 iterations to converge independent of the mesh size.

Concentration subproblem
Let us consider algebraic equation (8) and assume that a spacial configuration of the magnetic field H is given. The magnetic field is computed from a solution of the magnetostatic problem (13) and the integral condition (7) Here ω i is the area of a circular element, obtained by rotating of the cell T i . A system of M nonlinear equations (21) and one linear equation (22) F(x)=0, where F =(F 1 ,...,F M+1 ) T , x =(C 1 ,...,C M , c c ) T , will be solved by the Newton method Here Eliminating C k+1 1 ,...,C k+1 M from the system (23) we get Finally we compute C k+1 i from the i-th equation of the system (23) for the given c k+1 The computations for C k+1 1 ,...,C k+1 M can be realized in parallel. The Newton method requires 3 − 5 iterations to converge for the relative a-posteriori error |x k+1 − x k | to be smaller than ǫ c (generally 10 −9 ). We set C 0 i = C 0 at the initial Newton step. An initial value c 0 c is defined from the condition that far from the interface the concentration and the magnetic field intensity are known and equal C 0 and H 1 0 , respectively. From equation (20) we get

Free-surface subproblem
Let us write equation (12) for space variables dimensionless over a We assume that the magnetic field H and the concentration C are given at the interface Γ. The magnetic field is determined from the fluid side as H = |∇φ 1 | and H n = ∇φ 1 · n for the normal vector n, external to the fluid domain Ω 1 . We introduce dimensionless parameters

Equation (25) takes a new form
where and the effective field The unknown constant c f is different in equations (25) and (26) and will be fixed later. The integrand in (27) is further transformed We approximate the integral term in (27) by a composite trapezoidal rule on a uniform grid Here h = H/n, (n + 1) is the number of the grid points, H i = ih and C i denotes the concentration, corresponding to the field H i . f I (C 0 ,0)=0, because L(0)=0andL ′ (0)=1/3. The final form of the function f (C, H), used in computations, is To find concentration C i , corresponding to the integration points H i , we apply the Newton method to equation (20) for the given field H = H i and the given c c . The value of constant c c is defined in the process of solving the concentration subproblem. For details of computations see Section 3.2. Test computations of the integral approximation (28) for n = 10 and n = 20 show changes in the fifth significant digit. We use n = 20 for our computations.
The free boundary Γ is represented by an arc-length parametrization where the parameter s is measured from the top of the peak (s = 0) to the peak foot (s = ). Following the approach in (Polevikov, 2004), we reformulate equation (26) as a system of second-order ordinary differential equations Herer ( . The constant c f is determined by integrating equation (30) Equations (30) are augmented by boundary conditions The nonlocal boundary condition is due to the integration by parts of the volume conservation condition (31).
An iterative finite-difference scheme of the second order approximation for the parametric Young-Laplace equations was developed in (Polevikov, 2004). We apply the developed approach to equations (30) with boundary conditions (32) Here {r i } N i=0 and {z i } N i=0 are grid-functions uniformly distributed over the free surface with a step size h = 1/N. The difference quotients correspond to the central derivatives (r• s ,z• s )and the second derivatives (r ss ,z ss ). Nonlinearities of equations (30) are resolved by iterations, resulting in a three-diagonal system for the unknown grid functions at the (k + 1)-th iteration. A relaxation technique with a parameter τ f is applied to improve numerical stability. We took τ f = 0.01 in computations.

Decoupling strategy
The model under study consists of the magnetostatic subproblem, the concentration subproblem and the free-surface subproblem. The magnetostatic subproblem is described by a nonlinear elliptic partial differential equation with jumping coefficients (13) for the magnetostatic potential, augmented by transmission and Dirichlet-Neumann boundary conditions (14)-(16). The concentration subproblem is presented by nonlinear algebraic equation (20) for the concentration, augmented by integral condition (7). The free-surface subproblem is described by a system of two nonlinear ordinary differential equations (30) for the parametric representation of the free surface, augmented by integral and boundary conditions (31)-(32).
We apply an iterative decoupling strategy for solving the coupled system of equations. It consists of three steps at every iteration. The first step is a numerical solution of the magnetostatic problem in a fixed domain and for a given distribution of the concentration. The finite element method is applied for the discretization, see Section 3.1 for details. The second step is a numerical solution of system of nonlinear equations (21)-(22) for the element-wise concentration at the finite-element mesh for the given magnetic field distribution, as a solution at the first step. The Newton method is applied for a solution of the system, see Section 3.2 for details. The third step is a numerical solution of the free-surface subproblem for the given magnetic field from the first step and the given concentration from the second step of the iterative procedure. The finite-difference method is applied for the discretization, see Section 3.3 for details. A relaxation technique is applied to the free surface representation at every iteration It allows to suppress a rapid change of free surface shapes during iterations. We take initially τ = 0.1 and decrease this value to τ = 0.01 in the case of strong shape changes.
An initial surface configuration at the first iteration of the presented iterative algorithm is defined as a small perturbation of the flat surface with an amplitude of around 1 % of the cell radius. An initial concentration equals C 0 . The iterations are stopped when the change in the surface shape is smaller than a prescribed threshold ǫ ( generally 10 −7 ) The iterative process is controlled by the threshold ǫ, whereas three subproblems are controlled by own thresholds ǫ p , ǫ c and ǫ f .
All algorithms, discussed in Section 3, and the coupling of three subproblems were implemented in Fortran with the help of the software tools, earlier developed for the Rosensweig instability computations. Numerical results of the previous computations are published in (Bashtovoi et al., 2002;Lavrova et al., 2008;2010).

Results of computations
Numerical calculations were performed for the magnetic fluid EMG 901 (Ferrotec) with the following characteristic properties: the initial susceptibility χ = 2.2, the density ρ = 1406 kg/m 3 , the surface tension coefficient σ = 0.025 kg/s 2 , the saturation magnetization M s = 48 kA/m, the volumetric concentration, corresponding to the uniform particle distribution, C 0 = 0.1. The initial susceptibility of the Langevin magnetization χ L is related to the initial susceptibility of the effective-field magnetization (2) as see e.g. (Pshenichnikov et al., 1996). For the considered magnetic fluid we have that χ = 2.2 and χ L ≈ 1.47489. The control parameter of the model is the applied field intensity H 0 .
Computations are performed at the computational domain Ω ax =(0, a) × (−δ, δ) with δ = 5a. Computations with different δ have shown that the error caused by replacing the unbounded domain by a bounded one is less than 1%.
A linear stability analysis for the Rosensweig instability was carried out under the assumption of a uniform particles distribution in a layer of infinite thickness (Rosensweig, 1998). The stability theory predicts a critical value of the magnetic field intensity H c , corresponding to the onset of instability, as a solution of the nonlinear equation The intensity H c corresponds to the fluid domain. The critical field intensity at the air domain H * is found from the transmission condition [B · n]=0, satisfied at the unperturbed interface We get H * ≈ 9.11 kA/m for the considered ferrofluid. The stability theory predicts a critical value of the pattern wavelength λ c = 2π/ ρg/σ.
We assume that the hexagonal pattern wavelength λ hex , see Fig. 3, equals λ c .T h e nf o rt h e radius a of the circular cell we have We assume that the parameter λ is fixed for any applied field intensity.
Two meshes have been used for computations to analyze an influence of the discretization refinement to the computational predictions. Table 1 shows the critical field, the maximum value of the particle concentration over the fluid domain and z-coordinate of the equilibrium free surface at the peak axis and the peak foot for the applied field H 0 = 9.2 kA/m at different meshes. The found difference in values allows us to conclude that computations at the mesh with 160 × 1600 nodes are accurate enough. This mesh has been used to get results in Fig. 4-Fig. 6.
Mesh Results of two models, which account for a nonuniform particle distribution inside the ferrofluid layer will be compared with the results of the model with a uniform particle distribution. The first model assumes no interaction between particles and it was numerically studied in (Lavrova et al., 2010). The second model accounts for interaction between particles and is the subject of this contribution. The model with a uniform particle distribution is called model 3 in what follows.
Computations for the first and the third models in (Lavrova et al., 2010) show that the onset of the instability occurs at H * 1 = 9.12 ± 0.01 kA/m and H * 1 = H * 3 . This value nearly coincides with the result of the linear stability analysis H * ≈ 9.11 kA/m, which assumes a uniform particle distribution. It means, that the concentration effect does not influence to the onset of the instability in the frame of the model without particle interactions. Computations for the second model, which takes into account particle interactions, show that the onset of the instability occurs in a weaker field H * 2 = 8.65 ± 0.01 kA/m. A concentration effect in this case influences to the critical field. A possible reason for this effect is that the interparticle interaction can intensify considerably the fluid mangetization and a small initial surface perturbation is developed to a surface pattern in a weaker field if to compare with the models without particle interactions. Fig. 4 shows equilibrium free-surface shapes for three models. A more elongated peak region is formed for solutions with a nonuniform particle distribution. Namely, the peak is 20 % higher for the model without particle interactions and 60 % higher for the model with particle interactions in comparison with the uniform case. A 20 % difference of the first model is due to the concentration effect, which intensifies a spacial nonuniformity of the fluid magnetization in the peak region. A 60 % difference of the second model is influenced also by the fact that the onset of the instability for the second model occurs at the weaker field H * 2 < H * 3 .Af u r t h e r increase in the field strength results in the increasing peak amplitude, which can lead to a sizable difference in the peak height if we compare results of the second and the third model at the overcritical field H 0 = 9.2 kA/m. Fig. 4 contains isolines of the equilibrium concentration for the second model. The main inhomogeneity of the particle distribution occurs at the peak region. The concentration takes the greatest value at the top of the peak and the smallest value at the peak foot. Fig. 4. Overcritical equilibrium free-surface shapes at the applied field H 0 = 9.2 kA/m. 1 -nonuniform particle distribution without particle interaction, 2 -nonuniform particle distribution with particle interaction, 3 -uniform particle distribution. Isolines of the concentration corresponds to C/C 0 = {0.995, 0.999, 1.001, 1.02, 1.05, 1.08, 1.1, 1.11}. Fig. 5 shows the equilibrium distribution of the particle concentration over the peak axis for three models. The concentration increases monotonically in z-direction, moving along the peak axis, for the models with the nonuniform particle distribution and the concentration is constant for the third model. The concentration takes a value at the peak top which is about 25 % greater than in the fluid bulk for the first model and about 11 % greater for the second model. Taking into account the particle interaction, we get a smaller concentration at the peak but a more elongated shape. Fig. 5 shows that the concentration equals the volumetric value C 0 for z/a < −1 and the particle diffusion mechanism is present only near the free surface. Fig. 5. Equilibrium distribution of the particle concentration over the peak axis at the applied field H 0 = 9.2 kA/m. 1 -nonuniform particle distribution without particle interaction, 2nonuniform particle distribution with particle interaction, 3 -uniform particle distribution.
The distribution of a z-component of the magnetic field vector inside of the ferrofluid and in the air is presented in Fig. 6 for three models. The magnetic field is uniform inside of Fig. 6. Equilibrium distribution of the field intensity over the peak axis at the applied field H 0 = 9.2 kA/m. 1 -nonuniform particle distribution without particle interaction, 2nonuniform particle distribution with particle interaction, 3 -uniform particle distribution. the ferrofluid far from the interface, H z = H 1 0 , and the particle diffusion mechanism is absent there. The field distribution for three models coincides for z/a < −1. The intensity of the field increases monotonically in z-direction, moving along the peak axis inside the ferrofluid, up to the value close to the applied field intensity H 0 at the interface. Crossing the interface, the field intensity jumps to the value 3.17H 0 for the first model, 2.79H 0 for the second model and 2.7H 0 for the third model. The intensity of the field decreases monotonically in z-direction, moving along the peak axis inside the air, up to the value of the applied field intensity H 0 far from the interface. All considered models show that the field non-uniformity occurs at the region −a < z < 2a and the field is nearly uniform outside of this region. It means that a restriction of the unbounded domain by a bounded one with −5a < z < 5a will introduce an insignificant error in computations.

Conclusions
The effect of particle interaction in ferrofluid has been taken into account in numerical simulations of the Rosensweig instability for the first time. The mathematical model is based on a recently developed mass transfer equation for describing diffusion of interacting ferromagnetic particles in ferrofluids in (Pshenichnikov et al., 2011). The mathematical model for the Rosensweig instability in homogeneous ferrofluids is augmented by the concentration subproblem in the form of the nonlinear algebraic equation of the particle concentration and the magnetic field intensity inside the ferrofluid. We suggest an approach for numerical solution of the concentration subproblem and for the coupling of the concentration subproblem to the magnetic field and free surface computations.
Based on comparison between three models (model 1 -nonuniform particle distribution without particle interaction, model 2 -nonuniform particle distribution with particle interaction, model 3 -uniform particle distribution), it has been shown that the effect of particle interaction cannot be neglected or replaced by simpler models which try to capture the diffusion of particles. It was found that the onset of the instability occurs in a weaker field and the peak region is more elongated, if we take into account the particle interactions.
The effect of particle interaction on the pattern wavelength can not be numerically analyzed in the frame of the presented mathematical model. The pattern wavelength is fixed by the critical value obtained in the frame of the linear stability analysis. The experimental study of the aging variation of the Rosensweig instability in (Sudo et al., 2006) shows, however, that the pattern wavelength changes with time dramatically. Namely, the cell pattern gradually bifurcates at the constant magnetic field, and the number of spikes increases with time. That is why we plan for the future to consider a soliton-like surface configurations in the Rosensweig instability (Richter & Barashenkov, 2005). The pattern wavelength is a variable of this problem and the effect of particle interaction to the wavelength selection can be numerically analyzed. This book demonstrates applications and case studies performed by experts for professionals and students in the field of technology, engineering, materials, decision making management and other industries in which mathematical modelling plays a role. Each chapter discusses an example and these are ranging from wellknown standards to novelty applications. Models are developed and analysed in details, authors carefully consider the procedure for constructing a mathematical replacement of phenomenon under consideration. For most of the cases this leads to the partial differential equations, for the solution of which numerical methods are necessary to use. The term Model is mainly understood as an ensemble of equations which describe the variables and interrelations of a physical system or process. Developments in computer technology and related software have provided numerous tools of increasing power for specialists in mathematical modelling. One finds a variety of these used to obtain the numerical results of the book.