Next Article in Journal / Special Issue
Finite Element Analysis of Laminar Heat Transfer within an Axial-Flux Permanent Magnet Machine
Previous Article in Journal
Robust H Loop Shaping Controller Synthesis for SISO Systems by Complex Modular Functions
Previous Article in Special Issue
Quasi-Analytical Model of the Transient Behavior Pressure in an Oil Reservoir Made Up of Three Porous Media Considering the Fractional Time Derivative
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Riccardo Fazio
* and
Alessandra Jannelli
Department of Mathematical and Computer Science, Physical Sciences and Earth Sciences, University of Messina, Viale F. Stagno D’Alcontres 31, 98166 Messina, Italy
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2021, 26(1), 22; https://doi.org/10.3390/mca26010022
Submission received: 23 February 2021 / Accepted: 7 March 2021 / Published: 9 March 2021
(This article belongs to the Special Issue Advances in Computational Fluid Dynamics and Heat & Mass Transfer)

Abstract

:
This paper deals with a non-standard implicit finite difference scheme that is 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. The proposed approach allows imposing the given boundary conditions at infinity exactly. 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.

1. 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].

2. 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]
d 3 u d x 3 + u d 2 u d x 2 + β 1 d u d x = 0 u ( 0 ) = d u d x ( 0 ) = 0 , d u d x ( ) = 1 ,
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.

3. The Finite Difference Scheme

Without a loss of generality, we consider the class of BVPs
d u d x = f x , u , x [ 0 , ) , g u ( 0 ) , u ( ) = 0 ,
where u ( x ) is a d dimensional vector with u ( x ) for = 1 , , d as components, f : [ 0 , ) × I R d I R d , and g : I R d × I R d I R d . 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 quasi-uniform 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),
x = c · ln ( 1 ξ ) ,
and
x = c ξ 1 ξ ,
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 quasi-uniform 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 c N 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
c · ln ( 1 ξ ) < c ξ 1 ξ ,
for all ξ , but ξ = 0 .
The problem under consideration can be discretized by introducing a uniform grid ξ n of N + 1 nodes in 0 , 1 with ξ 0 = 0 and ξ n + 1 = ξ n + h with h = 1 / N , so that x n is a quasi-uniform grid in 0 , . The last interval in (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
x n + α = x ξ = n + α N ,
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
u n + 1 / 2 x n + 3 / 4 x n + 1 / 2 x n + 3 / 4 x n + 1 / 4 u n + x n + 1 / 2 x n + 1 / 4 x n + 3 / 4 x n + 1 / 4 u n + 1 , d u d x n + 1 / 2 u n + 1 u n 2 x n + 3 / 4 x n + 1 / 4 .
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
U n + 1 U n a n + 1 / 2 f x n + 1 / 2 , b n + 1 / 2 U n + 1 + c n + 1 / 2 U n = 0 , f o r n = 0 , 1 , , N 1 g U 0 , U N = 0 ,
where
a n + 1 / 2 = 2 x n + 3 / 4 x n + 1 / 4 , b n + 1 / 2 = x n + 1 / 2 x n + 1 / 4 x n + 3 / 4 x n + 1 / 4 , c n + 1 / 2 = x n + 3 / 4 x n + 1 / 2 x n + 3 / 4 x n + 1 / 4 ,
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
1 d ( N + 1 ) = 1 d n = 0 N | Δ U n | TOL ,
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.

4. 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 quasi-uniform 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
d 1 u d x = 2 u , d 2 u d x = 3 u , x ( 0 , ) d 3 u d x = 1 u 3 u β ( 1 2 u ) ,
with
1 u ( 0 ) = 2 u ( 0 ) = 0 , 2 u ( ) = 1 ,
or, in an equivalent form,
u = ( 1 u , 2 u , , 3 u ) T , f ( x , u ) = 2 u , 3 u , 1 u 3 u β ( 1 2 u ) T , g ( u ( 0 ) , u ( ) ) = ( 1 u ( 0 ) , 2 u ( 0 ) , 2 u ( ) 1 ) T ,
where u ( x ) is a three-dimensional vector with components u ( x ) for = 1 , 2 , 3 , and f : [ 0 , ) × I R d I R d and g : I R d × I R d I R d , with d = 3 .
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
1 u ( x ) = 0.5 x , 2 u ( x ) = 1 , 3 u ( x ) = exp ( x ) .
In Table 1, we report the numerical results that were obtained for the missing initial data d 2 u d x 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 d x 2 ( 0 ) = 1.177255155089504 , obtained in 6 iterations.
Table 2 and Table 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
R e l E r r = D T M F D F D R e l E r r = N I T F D F D
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 d x 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 = r N 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 d x 2 ( 0 ) . We can conclude that the reported extrapolated value is correct up to 9 decimal places.

5. 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 Table 2 and Table 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.

Author Contributions

The authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

The research of this work was supported, in part, by the University of Messina and by the GNCS of INDAM.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Blasius, H. Grenzschichten in Flüssigkeiten mit kleiner Reibung. Z. Math. Phys. 1908, 56, 1–37. [Google Scholar]
  2. De Hoog, F.R.; Weiss, R. An approximation theory for boundary value problems on infinite intervals. Computing 1980, 24, 227–239. [Google Scholar] [CrossRef]
  3. Lentini, M.; Keller, H.B. Boundary value problems on semi-infinite intervals and their numerical solutions. SIAM J. Numer. Anal. 1980, 17, 577–604. [Google Scholar] [CrossRef]
  4. Markowich, P.A. A theory for the approximation of solution of boundary value problems on infinite intervals. SIAM J. Math. Anal. 1982, 13, 484–513. [Google Scholar] [CrossRef]
  5. Markowich, P.A. Analysis of boundary value problems on infinite intervals. SIAM J. Math. Anal. 1983, 14, 11–37. [Google Scholar] [CrossRef] [Green Version]
  6. Kitzhofer, G.; Koch, O.; Lima, P.; Weinmüller, E. Efficient numerical solution of the density profile equation in hydrodynamics. J. Sci. Comput. 2007, 32, 411–424. [Google Scholar] [CrossRef]
  7. Fazio, R. The Blasius problem formulated as a free boundary value problem. Acta Mech. 1992, 95, 1–7. [Google Scholar] [CrossRef]
  8. Boyd, J.P. Chebyshev and Fourier Spectral Methods; Dover: New York, NY, USA, 2001. [Google Scholar]
  9. Shen, J.; Wang, L. Some recent advances on spectral methods for unbounded domains. Commun. Comput. Phys. 2009, 5, 195–241. [Google Scholar]
  10. Fazio, R.; Jannelli, A. Finite difference schemes on quasi-uniform grids for BVPs on infinite intervals. J. Comput. Appl. Math. 2014, 269, 14–23. [Google Scholar] [CrossRef]
  11. Fazio, R.; Jannelli, A. BVPs on infinite intervals: A test problem, a nonstandard finite difference scheme and a posteriori error estimator. Math. Methods Appl. Sci. 2017, 40, 6285–6294. [Google Scholar] [CrossRef]
  12. Fazio, R.; Jannelli, A.; Rotondo, T. Numerical study on gas flow through a micro-nano porous medium based on finite difference schemes on quasi-uniform grids. Int. J. Nonlinear Mech. 2018, 105, 186–191. [Google Scholar] [CrossRef]
  13. Fazio, R.; Jannelli, A. Two finite difference methods for a nonlinear BVP arising in physical oceanography. Atti Accad. Peloritana Pericolanti 2018, 96, A3. [Google Scholar]
  14. Fazio, R. Perpetual American put option: An error estimator for a non-standard finite difference scheme. arXiv 2014, arXiv:1412.1621. [Google Scholar]
  15. Singh, B.B.; Chandarki, I.M. Non-integral technique and differential transformation method for MHD boundary layer flow of an incompressible fluid past a flat plate. Int. J. Appl. Math. Res. 2012, 1, 46–64. [Google Scholar] [CrossRef] [Green Version]
  16. Brighi, B. The equation f”’ + ff” + g(f’) = 0 and the associated boundary value problems. Results Math. 2012, 61, 355–391. [Google Scholar] [CrossRef]
  17. Lambert, J.D. Numerical Methods for Ordinary Differential Systems; Wiley: Chichester, UK, 1991. [Google Scholar]
  18. Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A., Jr. Spectral Methods, Fundamentals in Single Domains; Springer: Berlin, Germany, 2006. [Google Scholar]
  19. van de Vooren, A.I.; Dijkstra, D. The Navier–Stokes solution for laminar flow past a semi-infinite flat plate. J. Eng. Math. 1970, 4, 9–27. [Google Scholar] [CrossRef]
Figure 1. Numerical solution for the problem (1) for β = 1.2 and N = 80 .
Figure 1. Numerical solution for the problem (1) for β = 1.2 and N = 80 .
Mca 26 00022 g001
Table 1. Numerical results that are related to d 2 u d x 2 ( 0 ) for different values of N.
Table 1. Numerical results that are related to d 2 u d x 2 ( 0 ) for different values of N.
NIter d 2 u dx 2 ( 0 )
106 1.1790527
206 1.1776846
406 1.1773411
806 1.1772552
1606 1.1772336
3207 1.1772282
6407 1.1772269
12807 1.1772261
Table 2. Numerical results, related to d 2 u d x 2 ( 0 ) , and comparison with results by [15].
Table 2. Numerical results, related to d 2 u d x 2 ( 0 ) , and comparison with results by [15].
β DTM [15]FD (This Study)Err Rel
0.00.469100.46958390.0010
0.20.663430.63899120.0382
0.40.800090.77497390.0324
0.60.916590.89175360.0279
0.81.019880.99563430.0244
1.01.113621.09008150.0216
1.21.200061.17724480.0194
1.41.280681.25856680.0176
1.61.356521.33507100.0161
1.81.428341.40751420.0148
2.01.496711.47647510.0137
Table 3. Numerical results, related to d 2 u d x 2 ( 0 ) , and comparison with results by [15].
Table 3. Numerical results, related to d 2 u d x 2 ( 0 ) , and comparison with results by [15].
β NIT [15]FD (This Study)Err Rel
0.00.469200.46958390.0008
0.20.648190.63899120.0174
0.40.787490.77497390.0162
0.60.905620.89175360.0156
0.81.010020.99563430.0145
1.01.104601.09008150.0133
1.21.191701.17724480.0123
1.41.272851.25856680.0113
1.61.349131.33507100.0105
1.81.421321.40751420.0098
2.01.490021.47647510.0092
Table 4. Extrapolated values at origin x = 0 for d 2 u d x 2 ( 0 ) with β = 1 .
Table 4. Extrapolated values at origin x = 0 for d 2 u d x 2 ( 0 ) with β = 1 .
N g 3 U g , 0 3 U g , 1 3 U g , 2
100 1.090081494 --
200 1.090069055 1.090064908 -
400 1.090065945 1.090064908 1.090064908
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fazio, R.; Jannelli, A. A Non-Standard Finite Difference Scheme for Magneto-Hydro Dynamics Boundary Layer Flows of an Incompressible Fluid Past a Flat Plate. Math. Comput. Appl. 2021, 26, 22. https://doi.org/10.3390/mca26010022

AMA Style

Fazio R, Jannelli A. A Non-Standard Finite Difference Scheme for Magneto-Hydro Dynamics Boundary Layer Flows of an Incompressible Fluid Past a Flat Plate. Mathematical and Computational Applications. 2021; 26(1):22. https://doi.org/10.3390/mca26010022

Chicago/Turabian Style

Fazio, Riccardo, and Alessandra Jannelli. 2021. "A Non-Standard Finite Difference Scheme for Magneto-Hydro Dynamics Boundary Layer Flows of an Incompressible Fluid Past a Flat Plate" Mathematical and Computational Applications 26, no. 1: 22. https://doi.org/10.3390/mca26010022

APA Style

Fazio, R., & Jannelli, A. (2021). A Non-Standard Finite Difference Scheme for Magneto-Hydro Dynamics Boundary Layer Flows of an Incompressible Fluid Past a Flat Plate. Mathematical and Computational Applications, 26(1), 22. https://doi.org/10.3390/mca26010022

Article Metrics

Back to TopTop