Received 30 January 2016; accepted 26 March 2016; published 29 March 2016
Fluid flow through irregular geometries (that is, geometries that do not conform to known coordinate lines) is encountered in many natural and industrial phenomena. Natural occurrences include flow in channels and in rivers over the curvilinear channel bottom, flow into a dam with a curvilinear free-surface, and groundwater and oil flow through the curvilinear earth layers and reservoirs. Industrial applications of flow through irregular domains include flow of air over wings and streamlined bodies, flow in curvilinear ducts, pipes and conduits, flow through turbines and propeller systems, and flow through artificial human and animal tissues,  . Solutions to, and analysis of initial and boundary value problems are important in the above-mentioned processes, and many involve solving the Navier-Stokes equations.
In solving boundary value problems that involve viscous fluid flow in two space dimensions, the Navier- Stokes equations are usually cast in vorticity-streamfunction form. The streamfunction possesses Dirichlet conditions if the solid boundary is assumed to be a streamline of the flow. However, no explicit expression exists for the vorticity on a solid boundary. This forces one to impose a Neumann (derivative) boundary condition on the vorticity. Neumann boundary condition is derived using the definition of the vorticity as the Curl of the velocity vector,    . Vorticity at the solid boundary therefore remains a quantity to be solved for. Accurate computation of boundary vorticity has a great influence on the local and global accuracies of the numerical solution of quantities in the flow domain. Furthermore, local accuracy of the computation of boundary vorticity affects the overall global accuracy of the numerical schemed used. This emphasizes the urgent need of accurately computing the vorticity at the boundary  .
In order to capture the boundary effects more accurately, the use of non-uniform grid with grid clustering at the boundary proved to be a viable tool. Grid refinement has been extensively studied and widely used, and a fairly large number of algebraic and mapping methods have been proposed to produce grid clustering for any given boundary shape,      . However, grid clustering also produces unfavourable effects on the convergence and the global accuracy of numerical solution obtained. These effects are due in part to the altering of the coefficient matrix that represents the finite difference form of the governing equations, and the potential loss of diagonal dominance of the coefficient matrix,   . Additionally, non-smoothness characteristics of non-uniform grid have been shown to adversely affect the performance of the numerical scheme employed  .
The above conclusion has been reached for the case of uniform grid. However, many authors have reached a similar conclusion when a non-uniform grid is used, and some,  , emphasized that the accuracy of the solution to the governing equations depends on the order of the scheme and on other factors such, as the stretching parameters employed, the error tolerance used in the solution methodology, and the iterative procedure employed. In addition, oscillations in the solutions have been pointed out by Pandit, Kalita and Dalal,  , who stated that “For fourth-order approximation of the vorticity at boundary points, some spurious oscillations were observed. When third-order approximation was used for the same, initially slight oscillations were seen.”
In summary, computation of boundary vorticity accurately is tantamount to finding a most suitable finite difference expression for the first derivative of the tangential velocity component; one that minimizes the unfavourable effects on the coefficient matrix and one that produces the best approximation to the boundary vorticity, whether one employs a uniform or a non-uniform grid.
In the case of fluid flow in curvilinear domains, the domain is first mapped onto a rectangular domain (or a regular domain) where finite differences can be applied. Numerical solution is then sought for transformed governing equations and boundary conditions,   -      . In particular, the vorticity Neumann boundary condition is transformed into the new coordinate system and a suitable differencing scheme is employed to accurately approximate the boundary vorticity. Numerical approximation to the boundary vorticity is influenced by a number of factors, including:
- The type of grid and spacing used.
- Entry condition (say to a channel).
- The use of flow-field grid points.
In order to shed some light on the above factors, Awartani, Ford and Hamdan,  , and Siyyam,Ford and Hamdan,  , considered the flow of a viscous fluid in a two-dimensional curvilinear channel. The curvilinear physical domain is mapped onto a rectangular computational domain using the von Mises coordinates,     , and the governing equations and boundary conditions are cast in von Mises variables. Vorticity condition on the physical boundary is expressed in terms of the first derivative of the square of the speed of the flow on the boundary of the transformed plane. Their use of the von Mises coordinates reduces the effects of non-smoothness in coordinate clustering. This is due to the fact that when a uniform grid is chosen in the physical domain, use of the von Mises transformation produces a naturally smooth clustered grid in the computational domain. In the above previous work,   , schemes of first-, second- and third-order local accuracy were developed for the first derivative and used in the calculation of vorticity on the flow domain boundary. In those schemes, stencils involving up to five computational grid points were used, and a comparison was offered between uniform and non-uniform grids. Evaluation of vorticity at a single corner point in the flow domain (where the exact solution is known) demonstrated the superiority of non-uniform grid when the von Mises transformation is used.
The above reported previous work,   , however, remained silent about the use of schemes of higher- order local accuracy, and the impact of using clustering in the physical domain on the computed vorticity on the boundary. This motivates the current work in which we develop and validate finite difference expressions, for the first derivative, that can be used in the accurate evaluation of vorticity on the flow domain boundary. Ultimately, finite difference expressions that yield the best approximation for the vorticity on the boundary are recommended for implementation in the full numerical solution of the channel flow problem described in this work. The intention is to develop a fourth-order scheme of local accuracy to approximate the vorticity on the boundary for both uniform and non-uniform grids using five grid-point stencils, and to address the following points when transforming a curvilinear physical domain into a rectangular computational domain:
- The type of grid spacing used: Typically, clustered grid near the boundary is expected to produce more accurate results. However, in some coordinate transformation procedures there are factors limiting how “compacted” the grid points can be. An important question to ask is whether a uniform grid in the computational domain can produce a clustered grid in the physical domain.
- Entry conditions to a channel: In case of a uniform entry condition to a channel, there is a correspondence between a grid in the computational domain and a grid in the physical domain, in the sense that clustering in the physical plane produces the same clustering in the computational domain. However, when a more realistic parabolic entry profile is used, the correspondence is lost and the use of clustered grid becomes necessary in the computational domain. We will consider the effects of a parabolic entry profile when the von Mises transformation is used with both uniform and non-uniform grids.
- Use of interior grid points in the approximation: A question that needs to be answered is how many internal grid points should be used in approximating the vorticity on the boundary, and how far from the boundary can grid points be taken for uniform and clustered grids. This work will demonstrate that a natural order of grid points that are closest to the boundary result in more accurate approximations to boundary vorticity.
In order to present our procedures and report the findings of this study, this work is organized as follows. In the following sections, we present a formulation of the problem of viscous fluid flow in a two-dimensional curvilinear configuration and transform the governing equations into von Mises variables. We present a discretization of the computational domain and derive third and fourth order accurate standard finite difference schemes for the first derivative and we test these schemes on a uniform grid and four clustered grids.
2. Problem Formulation
2.1. The Physical Domain
We consider a physical problem that gives rise to the need for finite difference expressions for the first derivative. Non-uniform grid arises in a number of ways, including clustering of coordinates and where coordinate transformation is involved. In the current work we consider the von Mises transformation which gives rise to the need for non-uniform grid.
Consider the flow in a two-dimensional, symmetric, dimensionless channel bounded below and above by solid, curvilinear boundaries. The channel is shown in Figure 1, and is described by: where is a known smooth function.
The flow of a viscous, incompressible fluid through the channel is governed by the equations of continuity and linear momentum, which take the following dimensionless forms, respectively, (cf.   ):
wherein, for two-dimensional flow, is the velocity vector field, p is the pressure, Re is the Reynolds number, , , and subscript notation denotes partial differentiation. Equations (1) and (2) represent a system of three scalar equations that are to be solved for the primitive variables u(x,y), v(x,y), and p(x,y), subject to the no-slip velocity condition on the solid lower and upper boundaries. At the inlet and exit of the channel, the normal velocity component vanishes, that is, while the tangential velocity component, u, can be taken as a function of y (assuming the channel is long enough to allow for this condition to be imposed at the exit).
The dimensionless continuity Equation (1) implies the existence of a dimensionless streamfunction, , such that:
Elimination of the pressure term from Equation (2) and the introduction of the dimensionless vorticity that is defined in terms of the velocity components as
Figure 1. Representative sketch: physical domain.
facilitate casting the flow equations in the following vorticity-streamfunction form:
In vorticity-streamfunction, the problem has been reduced from that of solving three equations in three primitive variables to that of solving two equations in the two unknowns and.
Appropriate boundary and inlet conditions are the no-slip on the upper and lower walls of the channel, and a parallel, parabolic inlet profile at x = a. We assume that the channel is of sufficient length so that inlet and exit conditions are the same, namely:
The absence of explicit vorticity boundary conditions necessitates imposing appropriate inlet, exit and solid boundary conditions. Assuming that Equation (5) is valid on all boundaries then, with, can be used on all boundaries of the configuration at hand (even at the exit of the channel, since the profile of u is assumed to be known there). The above conditions translate into the following conditions on (obtained from Equations (3)-(4), (8a) and (8b)), and (obtained from Equations (5), (8a) and (8b)) along the channel inlet and exit, x = a, b:
Along the lower and upper boundaries, respectively, we have:
In order to solve Equations (6) and (7), subject to conditions (9a)-(9d), numerically using finite differences, the flow domain is first transformed into a rectangular domain using the von Mises coordinates, and the governing equations and boundary conditions are cast in von Mises variables. The von Mises transformation and the transformed governing equations and boundary conditions are discussed in what follows.
2.2. The Computational Domain and Von Mises Transformation
Consider the curvilinear net in which the curves = constant represent the streamlines of the flow. Now consider the transformation, defined by
The Jacobian of transformation is given by
If, then the inverse transformation exists and the following first partial derivative operators in the two coordinate systems are obtained,    :
Second partial derivative operators can be obtained by applying operators (12) and (13) onto themselves.
Applying the above transformation to the vorticity-streamfunction Equations (6) and (7), respectively, yields,    :
The problem has thus been transformed from that of solving Equations (6) and (7) for and, subject to conditions (9a)-(9d), in the irregular (x, y) physical plane to that of solving Equations (14) and (15) for and in the rectangular computational plane, shown in Figure 2, where = constant represents the streamlines of the flow, subject to the following transformed conditions,   :
on the lower boundary (17a)
on the upper boundary (17b)
at both inlet and exit for (17c)
at inlet and exit, and on the lower and upper boundaries. (17d)
The square of the speed of the flow is given by:
and the velocity components are expressed as
Clearly, vorticity condition on the solid upper and lower channel boundary is given by (17d) in terms of the first derivative of the square of the speed of the flow. In the process of solving Equations (14) and (15) numerically, subject to conditions (17), we need to develop finite difference formulas for the first derivative in (17d).
Figure 2. The transformed computational domain.
3. Finite Difference Expressions for the First Derivative
3.1. Discretizing the Flow Domain
The rectangular computational domain of Figure 2 is discretized using either a uniform grid or a clustered grid, with vertical grid lines ranging from I = 1 at x = a to I = Imax at x = b, and horizontal grid lines ranging from j = 1
at the lower computational boundary to j = Jmax at the upper computational boundary.
It has been shown,  , that if a uniform grid is fitted over the computational domain, the corresponding gridlines in the physical domain would be clustered away from the boundary and closer to the center of the channel. However, in order to capture the effects of boundary conditions more accurately, it is required to have the clustering near the boundary.
Clustering of grid lines near the computational boundary when the von Mises variables are used can be accomplished in two ways:
1) Selecting a uniform grid in the physical domain and calculating the grid spacing in the computational domain results in a “natural” clustering near the computational boundary. For the given parabolic inlet velocity profile, and the indicated range of streamlines, we can choose y to vary over the interval [−1,1], and calculate the step size between two adjacent grid lines. Suppose that we wish to produce a uniform grid of step size
in the physical domain, a clustered grid is computed using for. The first
three variable grid spacings in the computational domain, corresponding to a uniform grid in the physical domain with, are:, , , where
for j = 1, 2, 3. Clearly, the computational grid lines are clustered near the computational boundary, even though the physical grid lines were equally spaced (uniform grid). This points out the fact that the use of the von Mises variables has the added advantage of naturally producing clustered grid near the computational boundaries without resorting to an algebraic method to cluster the grid.
2) Selecting a grid in the physical domain that is non-uniform and clustered near the physical boundary produces a computational grid that is also clustered near the computational boundary.
Physical grid clustering can be accomplished in various ways, including the use of elementary functions. For the current work, we employ the square root, the cubic root and the square functions to define the clustering in the physical y-direction. These are illustrated in Tables 1-5, wherein we provide a comparison with the case of uniform grid in the physical domain (Table 1) and the case of uniform grid in the computational domain (Table 5) with. It should be noted that five grid points are chosen in this work (one boundary point, j = 1, and four internal points, j = 2, 3, 4, 5) since the development of a fourth-order accurate scheme requires a minimum of five grid points.
3.2. Calculating Square of the Speed at Inlet
At the inlet to the channel, the computational velocity conditions are:, , for. The square of the speed at inlet to the computational domain is given by
Table 1. Grid 1. Uniform grid in the physical domain for, 0.99, 0.98, 0.97, 0.96.
Table 2. Grid 2. Non-uniform grid in the physical domain. for, 0.99, 0.98, 0.97, 0.96.
Table 3. Grid 3. Non-uniform grid in the physical domain. for, 0.99, 0.98, 0.97, 0.96.
Table 4. Grid 4. Non-uniform grid in the physical domain. for, 0.99, 0.98, 0.97, 0.96.
Table 5. Grid 5. Uniform grid in the computational domain with.
Using (20), we compute for each of the five types of grid considered, at the first five grid points. These are illustrated in Table 6.
Table 6. Square of the speed at inlet to computational domain for the five grids used.
3.3. Derivation of Finite Difference Schemes
The definition of vorticity is expressed at the lower boundary of the computational domain, in von Mises coordinates as the following derivative of the square of the speed of the flow:
In order to obtain finite difference expressions for, we first develop forward finite difference schemes for the first partial derivative on the Right-Hand-Side of Equation (21) using boundary and internal grid points. To accomplish this, we will use the Taylor’s series expansion approach, and discuss both standard and non- standard schemes that can be developed.
In order to derive forward differencing schemes of local accuracy n for the first derivative along the grid line (i,1) using grid points, Siyyam, Ford and Hamdan,  , assumed the schemes to take the form:
where E is the local truncation error. The weights, , are calculated using Taylor’s series expansion of
about point (i,1), for, namely
The coefficients, , appearing in (23) and (24) are defined by
For uniform grid, = constant for all j, and (25) takes the form
Using (23) in (22), and equating to zero the coefficients of the first n partial derivatives (including the zero’th derivative) of with respect to, leads to the following n + 1 conditions on,:
The order of the leading error term, E, will be the coefficient of is given by:
Based on the above description, Siyyam, Ford and Hamdan,  , observed the following:
1) The relationship between the order of the scheme, n, and the minimum number, k, of grid points used is that. This will render system (28) determinate consisting of n + 1 equations in as many unknowns. The resulting difference scheme is referred to as a standard scheme. If is the maximum number of grid points used, then.
2) The number of standard schemes that use k = n + 1 grid points when grid points are employed, is
given by the combination:.
3) If, the number of unknown coefficients will be greater than the number of equations, and the resulting system will have an infinite number of solutions for the coefficients of the differencing scheme. In this case, assigning arbitrary values to the coefficients results in a non-standard or arbitrary scheme. In developing non-standard schemes, grid points are skipped. Skipping grid points results in deterioration of the solution obtained by the scheme.
4) If, the number of unknowns will be less than the number of equations and the system will have no solution. This implies that in a first-order scheme, n = 1, and k = n + 1 = 2 is the minimum number of grid points used. Using less than 2 grid points results in no scheme. Similarly, there is no second-order scheme that uses less than 3 grid points, and no third-order schemes that uses less than 4 grid points, and so on.
The above observations are used to construct Table 7, below, that gives the number of possible schemes when grid points are employed.
3.4. Derivation of Standard Fourth-Order Accurate Scheme
We employ the procedure above to derive a standard fourth-order accurate scheme (that is, n = 4) using five grid points. We thus take. The scheme will thus employ one boundary point and four internal grid points, namely grid points j = 1, 2, 3, 4 and 5. Assuming the scheme to be of the form
Table 7. Number of possible standard schemes corresponding to.
coefficients are selected to satisfy the following:
The leading local truncation error has the following leading term that is the coefficient of:
Equations (32) translate into the following five equations in the five coefficients:
Solution to this system of 5 equations in 5 unknowns is given by
For the four types of non-uniform grids given in Tables 1-4, the values of are summarized in Table 8.
Using (40)-(43) in (35)-(39), coefficients take the following forms in terms of non-uniform step-sizes,:
For uniform grid with step size h, coefficients take the following forms:
and the fourth-order expression for the first derivative takes the form:
Table 8. Values of for the non-uniform grids used.
Values of the coefficients for both uniform and non-uniform grids are given in Table 9, below.
3.5. Standard Third-Order Accurate Schemes
In order to arrive at a third-order (n = 3) locally accurate forward differencing scheme for the first derivative, Awartani, Ford and Hamdan,  , used the above procedure and employed four grid points () in a scheme of the form given by (22). Their derivation produced one standard scheme that uses grid points j = 1, 2, 3, 4. Their scheme is reproduced here for completeness and we derive other third-order accurate standard schemes.
then the coefficients satisfy the following conditions:
Equations (52) represent a determinate system consisting of 4 equations in 4 unknowns that has the following unique solution:
Table 9. Values of coefficients for the five grids used; fourth order accurate scheme.
The leading error term, , is the coefficient of.
The four standard, third-order accurate schemes that use the four grid points j = 1, 2, 3 and 4, are as follows:
Solution (53), above, provides the following expression that was initially derived by Awartani, Ford and Hamdan,  :
For uniform grid, (57) yields the following scheme:
Values of the coefficients for both uniform and non-uniform grids are given in Table 10, below.
Other standard third-order accurate schemes obtained in this work employ grid points and are of
Table 10. Values of coefficients of the (1,2,3,4) scheme.
the form given by (52), with coefficients. These are as follows:
For uniform grid, (59) yields the following scheme:
Values of the coefficients for both uniform and non-uniform grids are given in Table 11, below.
Table 11. Values of coefficients of the (1,2,3,5) scheme.
For uniform grid, (62) yields the following scheme:
Values of the coefficients for both uniform and non-uniform grids are given in Table 12, below.
For uniform grid, (65) yields the following scheme:
Table 12. Values of coefficients of the (1,2,4,5) scheme.
Values of the coefficients for both uniform and non-uniform grids are given in Table 13, below.
3.6. Standard First- and Second-Order Accurate Schemes
Six standard, second-order accurate schemes have been developed by Awartani, Ford and Hamdan,  , using combinations of grid points on the boundary and the four grid points closest to the boundary. Henceforth, the notation -scheme denotes the scheme that employs points on the grid lines For the sake of comparison of results, we list below the natural order (1,2,3)-scheme that produced the most accurate results in  :
For uniform grid, (3.49) reduces to the following second-order accurate formula:
Combinations of standard first-order schemes that use the two grid points have been derived by Siyyam, Ford and Hamdan,  , and take the following general form for (that is, when grid points are used):
For the sake of comparison of results, we list below the natural order (1,2)-scheme that produced the most accurate results in  :
For uniform grid, (72) reduces to:
Table 13. Values of coefficients of the (1,3,4,5) scheme.
4. Results and Discussion
An expression for the vorticity on the lower computational channel boundary is given by. The lower, left-corner vorticity expression is given by, with exact value of this corner vorticity
, as obtained from the channel inlet condition where y = −1.
In order to approximate, we approximate the first derivative of the square of the speed, , using
the schemes developed in this work. These values are computed on the uniform and non-uniform grids given in Tables 1-5, using the square of the speed values of Table 6.
In the following Tables, we list the computed values of together with the Percentage Relative Error (PRE) defined as:
Comparison of the schemes that use the natural order of grid points in conjunction with the five grids used, is given in Table 14. This table demonstrates that the fourth order accurate scheme produces a corner vorticity that is closest to the exact value, with a Percentage Relative Error of 0.125% when Grid 3 is used. The fourth order accurate scheme consistently produces lower Percentage Relative Error for all of the grids employed.
Standard, first-order accurate scheme of Equation (72) produces the results in Table 14, below, which shows that using Grid 3 produces the most accurate computed values of vorticity and least Percentage Relative Errors. This may be ascribed to the fact that Grid 3 uses grid points that are closest to the boundary, hence better able to capture boundary effects. Grid 5 (uniform grid) produces the least accurate computed values of vorticity and largest Percentage Relative Error. Contrasted with clustered grid, the uniform Grid 5 has its largest step-size (grid spacing) closest to the boundary, which explains the large Percentage Relative Error associated with it.
Standard, second-order accurate scheme that uses the natural order of grid points produces results illustrated in Table 14 which, consistently, emphasizes the superiority of the (1,2,3) scheme when using Grid 3. The Percentage Relative Error associated with this combination is 0.146%. The four, third-order accurate schemes produce results reported in Table 15, and echo similar conclusions to the superiority of the natural-order scheme and the use of Grid 3.
1) The Percentage Relative Error decreases with increasing accuracy of the scheme employed, for both variable and uniform grids.
2) Grids with lowest Percentage Relative Errors for all standard schemes are, respectively: Grid 3, Grid 2, Grid 1, Gird 4, Grid 5.
3) Schemes that employ the natural order of grid points produce smaller Percentage Relative Errors, for all grids employed.
Table 14. Comparison of values of and pre using standard first, second, third and fourth-order accurate schemes of natural order of grid points.
4) Using clustered Grid 3 produces the lowest Percentage Relative Error, for all standard schemes used.
5. Generalization of the Schemes to Computational Flow Field Grid Points
Forward-difference schemes that approximate the derivative at a typical grid point (i,j) in the computational flow domain, using up-to 5 grid points, employ the gridlines j, j + 1, j + 2, j + 3, j + 4, shown in Figure 3.
The standard forward-difference uniform and non-uniform schemes that are based on the natural order of grid lines are listed below for the sake of completeness and ease of reference.
(j, j + 1)-Scheme: first order
For non-uniform grid:
For uniform grid:
(j, j + 1, j + 2) Scheme: second order
For non-uniform grid:
For uniform grid:
(j, j + 1, j + 2, j + 3) Scheme: third order
Table 15. Values of and PRE using standard third-order accurate schemes.
Figure 3. Grid lines j, j + 1, j + 2, j + 3, j + 4, used in forward-difference schemes.
For non-uniform grid:
For uniform grid:
(j, j + 1, j + 2, j + 3, j + 4) Scheme: fourth order
For non-uniform grid:
For uniform grid:
For uniform grid with step size h, coefficients take the following forms:
and the fourth-order expression for the first derivative takes the form:
The main theme of this work has been the development and testing of numerical schemes for the first derivative using uniform and non-uniform grids. The schemes have been applied to the computation of corner vorticity in the case of viscous fluid flow through a two-dimensional curvilinear channel that has been mapped onto a rectangular computational domain using the von Mises transformation.
The following has been accomplished in this work.
1) A fourth-order accurate scheme that uses 5 grid points has been developed and tested.
2) Third-order accurate schemes not previously reported have been developed and tested in this work.
Results have clearly demonstrated that:
1) Non-uniform grids produce more accurate values of corner vorticity than uniform grids.
2) The clustered grid that consistently produced the best approximation to corner vorticity is Grid 3, which is based on a stretching of physical coordinates using the cubic root function. This produced the smallest grid- spacing closest to the boundary in both the physical and computational domains.
3) The schemes that produced the best approximation to corner vorticity are the standard schemes which employ the natural order of grid points.
S. O. Alharbi, gratefully acknowledges support for this work through his PhD scholarship from Majmaah University, Kingdom of Saudi Arabia.
*Corresponding author. On leave from Majmaah University, Kingdom of Saudi Arabia.
 Hamdan, M.H. (1997) Natural Coordinate System Approach to Coupled n-Phase Fluid Flow in Curved Domains. Applied Mathematics and Computation, 85, 297-304.
 Gupta, M.M. (1983) A Survey of Some Second-Order Difference Schemes for the Steady State Convection-Diffusion Equation. International Journal for Numerical Methods in Fluids, 3, 319-331.
 Spotz, W.F. and Carey, G.F. (1998) Formulation and Experiments with High-Order Compact Schemes for Non-Uniform Grids. International Journal of Numerical Methods for Heat & Fluid Flow, 8, 288-303.
 Agarwal, R.K. (1981) A Third-Order-Accurate Upwind Scheme for Navier-Stokes Solutions in Three Dimensions. In: Computers in Flow Predictions and Fluid Dynamics Experiments, Proceedings of the Winter Annual Meeting, Washington DC, 15-20 November 1981, 73-82.
 Hamdan, M.H. (2009) Recent Developments in the von Mises Transformation and Its Applications in the Computational Sciences: Plenary Lecture. In: Perlovsky, L., et al., Eds., Mathematical Methods, System Theory and Control, MAMECTIS09, 180-189.
 Shukla, R.K., Shukla and Zhong, X. (2005) Derivation of High-Order Compact Finite Difference Schemes for Non-Uniform Grid Using Polynomial Interpolation. Journal of Computational Physics, 204, 404-429.
 Siyyam, H.I., Ford, R.A. and Hamdan, M.H. (2009) First-Order Accurate Finite Difference Schemes for Boundary Vorticity Approximations in Curvilinear Geometries. Applied Mathematics and Computation, 215, 2378-2387.
 Castillo, J.E., Hyman, J.M., Shashkov, J.M. and Steinberg, S. (1995) The Sensitivity and Accuracy of Fourth Order Finite-Difference Schemes on Non-Uniform Grids in One Dimension. Computers and Mathematics with Applications, 30, 41-55.
 Pandit, S.K., Kalita, J.C. and Dalal, D.C. (2008) A Fourth Order Accurate Compact Scheme for the Solution of Steady Navier-Stokes Equations on Non-Uniform Grids. Computers and Fluids, 37, 121-134.
 Visbal, M.R. and Gaitonde, D.V. (2002) On the Use of Higher-Order Finite Difference Schemes on Curvilinear and Deforming Meshes. Journal of Computational Physics, 181, 155-185.
 Awartani, M.M., Ford, R.A. and Hamdan, M.H. (2005) Computational Complexities and Streamfunction Coordinates. Applied Mathematics and Computation, 169, 758-777.