Mathematical Analysis and Numerical Simulation in Magnetic Recording

The head-tape interaction in magnetic recording is described in the literature by a coupled system of partial differential equations. In this paper we study the limit case of the system which reduces the problem to a second order nonlocal equation on a one-dimensional domain. We describe the numerical method of resolution of the problem, which is reformulated as an obstacle one to prevent head-tape contact. A finite element method and a duality algorithm handling Yosida approximation tools for maximal monotone operators are used in order to solve numerically the obstacle problem. Numerical simulations are introduced to describe some qualitative properties of the solution. Finally some conclusions are drawn.


Introduction
In magnetic storage systems, data is written onto magnetic tapes driven over a magnetic head which generates a magnetic field, which encodes the data into the tape. In similar way, the magnetic head reads the data stored on the tape. The motion of the tape introduces a film of air into the narrow space between the head and the tape. The profile of the head is designed to minimize the spacing between the head and the tape to effect efficient signal transfer. For that purpose, the data storage industry introduces trenches into the head to control the tape position and to avoid contact.
The unknowns of the problem are the position of the tape "u" and the pressure of the air "p". Compressible Reynolds equation models the pressure of the air while the position of the tape satisfies a beam equation. In the trenches, Reynolds equation is not effective and constant pressure is assumed. We consider the interval [x i , x i + L i ] (for i = 1, n), where Reynolds equation is satisfied and the trench "i" is enclosed in the interval [x i + L i , x i+1 ] (for i = 1 . . . n − 1), where p = 1 (after normalization). Notice that x 1 and x n + L n describes the beginning and the end of the head profile which contains n − 1 trenches (see Figure 1). We study the limit case where the Poiseuille flow is neglected (i.e. we do not consider the second order terms in Reynolds equation) and just the Couette flow describes the behavior of the air film between the head and the tape. The problem is described by a second order equation with nonlocal terms given by with Dirichlet boundary conditions u(0) = u(L) = 0, (1.2) where δ describes the profile of the head. In [10] the authors study the case where δ does not present any trench, under the concavity assumption δ (x) < 0. This assumption is very restrictive, not only mathematically, but also physically because magnetic heads do not generally satisfy the concavity condition. Indeed, in order to reduce the effect of air entrainment trenches are dug into the head (see [5]). In [11] the problem is studied for any regular function δ when compressible Reynolds equation models the air pressure inside the trenches. We assume throughout the paper that In this work we present a numerical resolution of the model and some numerical simulations in order to illustrate the behavior of the solutions.
The article is organized as follows: In Section 2 we detail the modelling and the derivation of the equation (1.1). The results on the existence of solutions are briefly described for reader s convenience in Section 3. More details on the existence of solutions can be found in Tello [12]. In Section 4 we describe the numerical method of resolution which is based on an obstacle problem approach and a duality algorithm handling Yosida approximation tools for maximal monotone operators. We present some numerical simulations to show the behavior of the solutions to the system. The proof of uniqueness of solutions is still an open problem, however the numerical simulations suggest that uniqueness of solution is expected. The article ends with a final section including a discussion and some conclusions.

Modelling
We assume that the pressure of the air "p" satisfies the modified compressible Reynolds equation. Then, after nondimensionalization (see [4]) we obtain with Dirichlet boundary conditions The tape deflection u is given by a fourth order linear equation of Euler-Bernoulli Beam equation type (see [4], [9,Chap. 6] and [1]). The position of the tape "u" is given by the beam equation, (see [4], [9, Chap.6] for details).
The coefficients µ and k are defined by where E is the Young modulus, I is the inertia moment of the tape, T is the tension of the tape, ρ is the tape's density, V is the velocity of the tape, p a is the ambient pressure. The coefficient k describes the influence in the tape deflection of the external forces applied over the tape, the "pressure" in this case.
Dimensional analysis of equations are frequently used in equations describing fluid mechanics to simplify the systems. Typical magnitudes of the parameters are (see [10] and [9, Chap.6]): Note that the tape is moved over the head at high speed V , which makes the parameter 1. After nondimensionalization, the orders of magnitude of p, h and the size of the head L n − x 1 denoted by P , H and X are 1, a comparative argument of the magnitudes of the terms in (2.1) shows that From the physical point of view it means that the convection dominates diffusion in the equation, so the variation of the pressure caused by the Couette flow (described by the second order terms) is small compared with the influence of the Poiseuille flow (described by the first order term). As a consequence, the air pressure behavior is described by the first order term. Hence, In the same way we analyze the tape deflection: since U ∼ 1, P ∼ 1, L ∼ 10, µ ∼ 10 −3 and k ∼ 10 3 the analysis of the order of magnitude of the terms in (2.3) gives The inequality (2.6) shows that the influence in the tape deflection caused by tape's tension dominates the equation, i.e. the axial force responsible of tension and velocity of the tape dominates the system, and the influence of the curvature of the tape is neglected. Then, thanks to (2.6), equation (2.3) is simplified to a second order equation Notice that some of the boundary conditions in (2.2), (2.4) need to be dropped. We take From (2.5) we see that ph = const.
Since the second derivative of the solution is zero in the intervals (0, x 1 ), (x i +L 1 , x i+1 ) and (x n +L n , L), the equation (1.1) is equivalent to the following system of n equations with the appropriate boundary conditions: (2.9)

Existence of Solution
In this section we briefly describe the results concerning the existence of solution. The complete proof of the theorem can be found in [12].
For reader's convenience we present a sketch of the proof. We first extend δ to the interval [0, L] by a continuous and regular functionδ such thatδ → δ in ( We consider λ := (λ 1 , . . . , λ n ) ∈ I and the problem (3.1) A sub-and supersolutions method gives existence of solutions. A fixed point argument is used to obtain the existence of a unique solution u (λ, x) to the problem (3.1). We construct a function f : I → R n defined as follows We notice that f is a continuous and locally Lipschitz function. The most technical part of the proof is to prove the existence of a fixed point of f in I. The proof is based on a Schauder's fixed point theorem. Due to the nature of the approach, uniqueness of solution can not be obtained by this method. The proof ends taking limits when → 0.
In this section we shall describe the numerical method employed to solve numerically the problem (2.8)-(2.9). In order to compute the solution we reformulate the problem in terms of an obstacle one relative to the unknown h = u − δ, so that we can guarantee that h > 0, that is to say, the tape keeps above the head. We shall follow the ideas of [1,2], and apply a finite element method together with a duality algorithm handling Yosida approximation tools for maximal monotone operators. These techniques have already been successfully used for example in [7] and [8]. To be precise, we are going to deal with the following complementarity formulation of the problem (2.8)-(2.9) in terms of h: where −α is the second derivative of δ, in other words, α(x) = −δ (x). The above complementarity problem (4.1) can be equivalently formulated, in order to apply the duality method (see [3]), in terms of the indicator function (see [6] for details) of the convex set of admissible solutions I K = {h ∈ H 1 (Ω), h ≥ 0}, where Ω = (x 1 , x n + L n ). Then, the equivalent variational formulation to (2.8)-(2.9) in terms of the new variable h and the indicator function I I K of the convex set I K is as follows: Next, by applying some results relative to the subdifferential calculus to the convex function I I K (see [13]), we can write the above formulation in the form: The method proposed in [3] introduces a new unknown q(x) defined in terms of a positive parameter ω by Therefore, the problem (4.2) admits this new formulation: ∀ψ ∈ H 1 (Ω) and q ∈ ∂I I K (h) − ωh. As ∂I I K is a maximal monotone operator, the following equivalence holds (see [3] ω λ is the Yosida approximation of the operator (∂I I K − ωI) with parameter λ > 0.
In order to discretize the nonlinear problem (4.3) with respect to the coordinate x, we employ Lagrange linear finite elements. Let us define where E denotes a standard finite element interval and T ∆x a uniform grid with step ∆x (we follow a similar notation to the one used in [8]).
So, the discretized problem is formulated as follows: The above nonlinear problem is solved with the following iterative scheme: 1. First we initialize h ∆x,0 .
Mathematical Analysis and Numerical Simulation in Magnetic Recording 341 where the multiplier updating is carried out through the identity: In the implementation to solve the discretized problem we use a relaxation parameter ν ∈ (0, 1], so, in fact we have that The introduction of a relaxation parameter facilitates the convergence process of the iterative scheme (see [8] where this kind of technique is also employed). The convergence of the duality method is guaranteed for values such that ωλ = 0.5 (see [3]). In this case, we have that The computations of (∂I I K ) ω 1 2ω (r) can be obtained in a similar fashion than the one used in [8].

Numerical simulations
In this section we present some numerical results obtained with two different types of δ and for each δ we shall consider the case with no trenches, the one with one trench and finally with two trenches. We have considered k = 10 3 , L = [0, 9] x 1 = 4 and x n + L n = 5. In the one trench case, the trench is located from x 1 +L 1 = 4.7 to x 2 = 4.9, and the two trenches in the last case are located from x 1 + L 1 = 4.4 to x 2 = 4.55, and from x 2 + L 2 = 4.8 to x 3 = 4.9.
Note that the function considered as δ should be in accordance with the physics of the problem, so we are interested in considering functions with a concave profile (see Figure 1). The functions which we have considered in the simulations presented here are the following: Case 1: δ(x) = x(9 − x) − 0.5. Notice that in [0, 4) and in (5,9] the solution is a straight line as the second derivative is null there, and that u(0) = u(9) = 0. The reason why we have introduced the term -0.5 is because adding −0.5 we get that δ(0) = δ(9) = −0.5 < u(0) = u(9) = 0. The results are shown in Figure 2.
This kind of profile is considered in [1]. Note that with this choice, conditions (1.4)-(1.5) are satisfied. Note also that these conditions are only assumptions made to develop the mathematical analysis of the model carried out in [12]. The results obtained in this case are shown in Figure 3.
Note that in the two cases, the concavity condition δ < 0 is satisfied, but while in the first case, δ is constant, in the second one, δ depends on x.
In the simulations presented here, we have considered a tolerance for the error in the duality algorithm equal to 10 −4 , the parameter ν is chosen to be ν = 0.5 which offers a good performance, and q ∆x,0 is chosen to be identically zero. The step considered for the space discretization is ∆x = 5 · 10 −4 . The convergence of the duality algorithm is guaranteed because of the choice ωλ = 0.5 (see [3]). The choice of the parameter ω, in absence of an analytical solution, has been carried out in order to obtain meaningful physical solutions with a behavior as expected (see [1] for more details). Note that optimal choice of the parameters is a function of the solution of the problem. Feasible solutions are obtained with ω of order 10 3 . The number of iterations is an increasing function of ω. All the results presented in this study, the ones in Figures 2, 4 and 3 , have been obtained with ω = 10 3 . For a tolerance equal to 10 −4 , the number of iterations in all the simulations presented here is of the order of 3000. It is also observed that the number of iterations slightly decreases when increasing the number of trenches. In the simulations whose results we show in Figure 2, the iterations for no trenches are 3383, with one trench 3093 and with two trenches 2864. The same behavior is observed in the simulations carried out to obtain the results that are shown in Figure 3 (the number of iterations in this case are 3381, 3092 and 2864). This behavior could be due to the fact that in the gaps the formulation is simpler.
In Figure 2, we show the results obtained for the case in which δ (plotted on the left) has no trench, one trench and two trenches. Notice that the results obtained for u (on the right) in the case of no trench reach values of the order of 2260, when one trench is considered, then u reaches a maximum of the order of 1980, while in the case of two trenches the maximum value is of the order of 1800. One can also observe that in the case of two trenches, the values for u are smaller but do not decrease as fast at the end of the domain as do the ones obtained for the case of one trench or no trench.
In Figure 4 we show some simulations corresponding to δ with no trenches. On the left, we show δ (solid lines on the left top corner and on the left bottom corner) and the initial data u 0 = h 0 + δ considered to initialize the scheme of resolution in two different simulations (dotted lines on the left top corner and on the left bottom corner). On the right, we show the results for u corresponding to the two different initial data u 0 presented here, to which the method converges. One can observe that in both cases the result is the same. A large number of simulations have been made and we have always obtained uniqueness of solution for the same δ. In view of these numerical results, one could infer that uniqueness of solution to the problem studied here is more than possible.
In Figure 3 we show the results obtained with δ = 1+ 1 − (x − 4.5) 2 . From these simulations one can derive the same conclusions as the ones commented regarding Figure 2. Note that in order to obtain the numerical results shown in Figure 3, we have considered the same values for the parameters and the same location of the trenches as the ones used to obtain the results presented in Figure 2. In doing so we are able to compare properly the results.

Conclusions
In this article we consider a mathematical model to describe the tape deflection when it is driven over a magnetic head. We study the case where a given geometry of the head presents a finite number of trenches. In these regions, Reynolds equation is replaced by a constant pressure to model the behavior of the air. The modelling follows [9] and [10], where the problem is studied for a concave and regular head profile respectively. The paper focuses on the limit case = µ = 0, where the second order terms in the Reynolds equation are neglected and the phenomenon is described only by the transport terms. The fourth order beam equation is simplified to a second order equation.
A sketch of the proof of the existence of solutions is also presented. The proof follows a fixed point argument based on a sub-and super-solution method. In order to obtain appropriate sub-solutions we should impose some extra assumption in the geometry of the head: (1.4)-(1.7). These assumptions guarantee the existence of solutions. We also notice that (1.4) is a necessary assumption to obtain existence of solutions, while the rest of assumptions may be relaxed. See Tello [12] for details.
The numerical resolution of the problem follows an obstacle problem argument introduced to avoid head-tape contact. From the numerical simulations one could derive as a conclusion that uniqueness of solutions to the problem may be possible (see Figure 4), nevertheless an analytical proof of uniqueness remains open. Numerical simulations show that the distance between the tape and the head decreases when trenches are introduced in the head profile δ. At the end of the head, according to the numerical results, the head-tape distance decreases at a lower rate when the number of trenches is increased. So it seems that more trenches, mainly at the end, guarantee that the tape keeps a proper distance to the head.