A Non-Standard Finite Difference Scheme for Magneto-Hydro Dynamics Boundary Layer Flows of an Incompressible Fluid Past a Flat Plate

This paper deals with a non-standard finite difference scheme defined on a quasi-uniform mesh for approximate solutions of the Magneto-Hydro Dynamics (MHD) boundary layer flow of an incompressible fluid past a flat plate for a wide range of the magnetic parameter. We show how to improve the obtained numerical results via a mesh refinement and a Richardson extrapolation. The obtained numerical results are favourably compared with those available in the literature.


Introduction
The simplest example of the application of the boundary layer theory is related to the celebrated Blasius [1] problem. This problem describes the flow around a very thin flat plate.
The first goal of this paper is to numerically solve, with great accuracy, the Magneto-Hydro Dynamics (MHD) boundary layer equation that governs the flow of an incompressible fluid past a flat plate by an implicit non-standard finite difference scheme that is defined on a quasi-uniform mesh that allows imposing the given boundary conditions exactly.
Numerical methods for problems, like the one considered in this paper, can be classified according to the numerical treatment of the boundary condition imposed at infinity. The oldest and simplest treatment is to replace infinity with a suitable finite value, the so-called truncated boundary. However, being the simplest approach, this has revealed, within the decades, some drawbacks that suggest not to apply it, especially if we have to face a given problem without any clue on its solution behaviour. Several other treatments have been proposed in the literature to overcome the shortcomings of the truncated boundary approach. In this research area, some treatments are worthy of consideration: the formulation of so-called asymptotic boundary conditions by de Hoog and Weiss [2], Lentini and Keller [3], and Markowich [4,5]; the reformulation of the given problem in a bounded domain, as first studied by de Hoog and Weiss and more recently developed by Kitzhofer et al. [6]; the free boundary formulation that was proposed by Fazio [7], where the unknown free boundary can be identified with a truncated boundary; the treatment on the original domain via pseudo-spectral collocation methods, see the book by Boyd [8] or the review by Shen and Wang [9] for more details on this topic; and, finally, a non-standard finite difference scheme on a quasi-uniform grid that is defined on the original domain by Fazio and Jannelli [10]. The proposed non-standard finite difference scheme has been successively modified by Fazio and Jannelli [11] and used in [12,13]. Moreover, this method has been applied to the numerical solution of the perpetual American put option problem of financial markets, see Fazio [14].
The key advantage of the presented approach is that the given BVP is solved on a semi-infinite interval by introducing a stencil that is constructed in such a way that the boundary conditions at infinity are exactly assigned. In this way, the whole infinite domain is taken into account in the mapping where the grid-points are located at the mid-point of each sub-interval, and then the difficulty that is caused by numerical treatment of the last infinite sub-interval is avoided.
Finally, via a mesh refinement, by using Richardson's extrapolation, we improve the accuracy of the obtained numerical results, showing that they are in excellent agreement with the solutions found in the literature. This study concludes by comparing the current numerical results with those that are given by the integral approximation method (ITM) and the non-integral technique (NIT) used by Singh and Chandarki [15].

Model Problem
We consider a steady two-dimensional flow of a viscous fluid on a flat plate in the presence of a given transverse magnetic field with small electric conductivity and a large transverse magnetic field. When introducing appropriate similarity variables, the governing equations can be reduced to the following boundary value problem (BVP) [15] where β is the magnetic parameter. Let us notice that, for β = 0, the BVP (1) reduces to the celebrated Blasius problem [1]. For a general class of problems, including (1), Brighi [16] obtained results on the existence, uniqueness, and boundedness of solutions.

The Finite Difference Scheme
Without a loss of generality, we consider the class of BVPs Here, and in the following, we use Lambert's notation for the vector components [17, pp. 1-5].
In order to solve the problem (2) on the original domain, we first discuss quasiuniform grids maps from a reference finite domain and then introduce, on the original domain, a non-standard finite difference scheme that allows for us to impose the given boundary conditions exactly. Let us consider the smooth strict monotone quasi-uniform maps x = x(ξ), the so-called grid generating functions, see Boyd [8, pp. 325-326] or Canuto et al. [18, p. 96 and where ξ ∈ [0, 1], x ∈ [0, ∞], and c > 0 is a control parameter. Accordingly, a family of uniform grids ξ n = n/N defined on interval [0, 1] generates one parameter family of quasiuniform grids x n = x(ξ n ) on the interval [0, ∞]. The two maps (3) and (4) are referred to as logarithmic and algebraic maps, respectively. As far as the authors' knowledge is concerned, van de Vooren and Dijkstra [19] were the first to use these kinds of maps. We notice that more than half of the intervals are in the domain with a length approximately equal to c and x N−1 = c ln N for (3), while x N−1 ≈ cN for (4). For both maps, the equivalent mesh in x is nonuniform, with the most rapid variation occurring with c x. The logarithmic map (3) gives slightly better resolution near x = 0 than the algebraic map (4), while the algebraic map gives much better resolution than the logarithmic map as x → ∞. In fact, it is easily verified that The problem under consideration can be discretized by introducing a uniform grid (3) and (4), namely [x N−1 , x N ], is infinite, but the point x N−1/2 is finite, because the non integer nodes are defined by with n ∈ {0, 1, . . . , N − 1} and 0 < α < 1. These maps allow for us to describe the infinite domain by a finite number of intervals. The last node of such grid is placed on infinity, so right boundary conditions are correctly taken into account.
We approximate the values of the scalar variable u(x) and its derivative at mid-points of the grid x n+1/2 , for n = 0, · · · , N − 1, while using non-standard difference discretizations We emphasize that the key advantage of our non-standard finite difference formulation is to overcome the difficulty of the numerical treatment of the boundary conditions at infinity. In fact, the formulae (5) use the value u N = u(∞), but not x N = ∞ and, then, the boundary conditions at infinity are taken into account in a natural way. For the class of BVPs (2), a non-standard finite difference scheme on a quasi-uniform grid can be defined using the approximations given by (5) above, and it can be written, as follows for n = 0, 1, . . . , N − 1. The finite difference formulation (6) has an order of accuracy O(N −2 ). It is evident that (6) is a nonlinear system of d (N + 1) equations in the d (N + 1) unknowns U = (U 0 , U 1 , . . . , U N ) T . For the solution of (6), we can apply the classical Newton's method along with the simple termination criterion where ∆ U n , for n = 0, 1, . . . , N and = 1, 2, . . . , d, is the difference between two successive iterate components and TOL is a fixed tolerance.

Numerical Results and Comparison
In this Section, we present the numerical results that were obtained by solving the mathematical model (1) using the non-standard finite difference scheme (6) on the quasiuniform grid that was defined by the logarithmic map (3) with control parameter c = 2. The control parameter c defines the grid points distribution in the original physical infinite domain and it is chosen in such a way that [x 0 , x N−1 ] is not too large, being x N → ∞. The chosen value c = 2 leads to obtaining a finer grid points distribution close to x = 0 and then a better resolution where the solution has a transient behavior, with an increasing spatial resolution going toward infinity. Now, let us rewrite the model (1) as a first-order system, as follows where u(x) is a three-dimensional vector with components u(x) for = 1, 2, 3, and For all tests, we consider the control parameter c = 2, a fixed tolerance TOL = 10 −8 , and, as a first guess for Newton's iteration, the following initial data In Table 1, we report the numerical results that were obtained for the missing initial data d 2 u dx 2 for increasing points number N and β = 1.2. Here, "iter "stands for the number of the Newton's iterations. It is well known that Newton's method could have a computational cost that is quite high, depending on the initial iterate and the analytic form of the function whose root is being sought. A suitable choice of the initial iterate can contribute to reducing the computational cost. In this context, the chosen initial guess leads to obtain a low iterations number: 6 or 7 iterations, depending on the values of the grid point number. In Figure 1, we show the numerical solution that was obtained for β = 1.2 and N = 80. The recovered value of the second order derivative of the solution at the origin is d 2 u dx 2 (0) = 1.177255155089504, obtained in 6 iterations.  Tables 2 and 3 list the obtained numerical results for different values of parameter β, with β = 0, 0.2, · · · , 2 and N = 100. For the sake of brevity, we have chosen to only report the values of the wall shear stress, which is the second derivative value at the origin. Within the same table, we can compare our results with those that were reported by Singh and Chandarki [15]. The last column, named Err Rel, reports the values of the relative errors, as computed by where, DTM, NIT, and FD represent the numerical solutions obtained, for different values of β, by the integral approximation method (ITM) and the non-integral technique (NIT) used by Singh and Chandarki [15], and by the proposed numerical method, respectively. In both cases, the value of the relative errors decrease as the values of β increase. Note that the relative errors are all of same the order and, moreover, for β = 0, we obtain the smaller relative error. This is due to the fact that β affects the non-linearity of the model under study. The problem with β = 0 corresponding to the Blasius problem and, in this case, the computed missing initial condition d 2 u dx 2 (0) = 0.4695839 can be compared with the value 0.469599988361 that wascomputed by Fazio [7] by a free boundary formulation of the Blasius problem. We apply Richardson's extrapolation, using several refinements of the computational domain, in order to improve the accuracy of the computed solution. On the computational domain of the problem, we build a quasi-uniform grid with a mesh-points number equal to N 0 and proceed with subsequent grid refinements by constructing meshes with grid-point numbers N g for g = 1, 2, · · · , where N g+1 = rN g with refinement factor r = 2. On each grid, the numerical solution U g , g = 0, 1, · · · , G is computed using the non-standard finite difference method. We adopt a continuation strategy in order to reduce the calculations; in fact, we use the final solution U g that was obtained on the grid g as an initial guess for calculating the solution U g+1 on the grid g + 1, where the new grid values are approximated by linear interpolations. We define the level of the Richardson's extrapolation by the index k and the two numerical solutions related to the grids g and g + 1 at the extrapolated level k by U g,k and U g+1,k . We use the following formula to calculate a more accurate approximation U g+1,k+1 = U g+1,k + U g+1,k − U g,k 2 p k − 1 k = 0, 1, · · · , G − 1 .
In Table 4, we report the extrapolated values with N = 100, 200, 400 grid points for β = 1. The last extrapolated value is 3 U 2,2 = 1.090064908 and it can be considered as our benchmark value for d 2 u dx 2 (0). We can conclude that the reported extrapolated value is correct up to 9 decimal places. Table 4. Extrapolated values at origin x = 0 for d 2 u dx 2 (0) with β = 1.

Concluding Remarks
In this paper, the mathematical model, which describes the MHD boundary layer flow of an incompressible fluid past a flat plate, is solved by the non-standard implicit finite difference method on a quasi-uniform grid for the different magnetic parameters β. We solved the given BVP on a semi-infinite interval by introducing a stencil that is built in such a way that the boundary conditions at infinity are exactly assigned. Moreover, we improved the accuracy of the numerical results by Richardson's extrapolation, which allowed for us to find a more accurate solution.
The values of the second-order derivative of the solution at the origin, for different values of parameter β, are reported in the Tables 2 and 3. The results are compared with those by Singh and Chandarki [15] in order to verify the accuracy of the proposed method. Moreover, in the case of the Blasius problem, we compare the missing initial condition with the one that was computed by Fazio [7] by a free boundary formulation of the Blasius problem. The computed values are found to be really accurate.