A CHEBYSHEV SPECTRAL METHOD FOR HEAT AND MASS TRANSFER IN MHD NANOFLUID FLOW WITH SPACE FRACTIONAL CONSTITUTIVE MODEL

In some recent studies, it has been suggested that non–Newtonian fluid flow can be modeled by a spatially non–local velocity, whose dynamics are described by a fractional derivative. In this study, we use the space fractional constitutive relation to model heat and mass transfer in a nanofluid. We present a numerically accurate algorithm for approximating solutions of the system of fractional ordinary differential equations describing the nanofluid flow. We present numerically stable differentiation matrices for both integer and fractional order derivatives defined by the one–sided Caputo derivative. The differentiation matrices are based on the series expansion of the unknown functions using a truncated Chebyshev polynomial of the first kind and interpolation using Gauss–Lobatto quadrature. We show that the proposed technique is highly effective for solving the fractional model equations.


INTRODUCTION
The boundary layer flow of both Newtonian and non-Newtonian fluids has been explored extensively in the literature because of the many industrial and scientific uses of these fluids. Experimental and theoretical investigations of non-Newtonian fluids have shown that they exhibit anomalous behaviours and are generally more complex than viscous fluids. To describe these behaviours, several models have been proposed to mathematically describe and experimentally quantify the stress-strain relationship. These constitutive relations include but are not limited to the power-law, Ellis and Powell-Eyring fluids models. These models and several others have been subjected to experimental and theoretical investigations. In some recent studies, attempts have been made to reconstruct the constitutive relations to incorporate hereditary properties, memory retention and non-local properties of non-Newtonian fluids. The primary objective is to recast the velocity gradient in the momentum transport equation as a fractional derivative.
Fractional calculus can best describe physical models with long term memory and non-locality of the spatial domain. For instance, constitutive relations with time fractional derivative are used to describe memory dependence in dynamical systems (Tarasov, 2013). An investigation and summary of studies on fluid flow problems in which the rheological properties of the fluid are described using fractional calculus is given by Sun et al. (2018). In this study, the classical constitutive equation of a viscous fluid was reconstructed so that the velocity gradient is expressed using a non-local spatial derivative. This approach allows for the capture of † Corresponding author. Email: shina@aims.edu.gh the non-locality of non-Newtonian fluid and the correlation between the fluid particles. They conducted a theoretical analysis of a steady pipe flow of a laminar incompressible flow, obtaining the frictional head loss, velocity profile and yield stress in terms of the fractional Reynolds number. Chen et al. (2015) investigated the fractional boundary layer flow of a Maxwell fluid on an unsteady stretching surface. The constitutive equation was recast so as to introduce space-time dependent fractional derivatives. Research on classical boundary layer flow, heat and mass transfer in magnetohyrodynamics nanofluid has been carried out by, among others, Dhanai et al. (2015) and Khan and Pop (2010). By applying the implicit finite difference method, velocity, temperature and nanoparticle concentration profiles were obtained in the study of Dhanai et al. (2015) which investigated the combined effects of heat and mass transfer parameters on a steady boundary layer nanofluid flow and examined multiple solutions of the conservation equations. Critical review of literature shows that heat and mass transfer in viscous nanofluid with fractional constitutive relation has not been studied. In this study, a fractional system of equations describing heat and mass transport in an incompressible magnetohydrodynamic viscous nanofluid is developed and solved using a spectral technique that uses Chebyshev polynomials of the first kind as basis functions.

MATHEMATICAL FORMULATION
We consider a steady two-dimension boundary layer flow of an incompressible nanofluid with space fractional shear rate induced by an im-permeable horizontal wall along the plane y = 0, where the y-axis is perpendicular to the wall and the x-axis along the wall. It is assumed that the flow occurs at y ≥ 0 and a constant magnetic field B is applied in the y-direction. The temperature and nanoparticle volume fraction at the wall (Tw, Cw) are constant and the ambient temperature and concentration are T∞ and C∞ respectively as in Fig. 1. By defining the shear stress associated with a fluid with fractional dynamic viscosity coefficient µα (kg/m 2−α /s) in the following form where α is the order of the fractional space derivative, and ∂ α u/∂y α denotes the fractional velocity gradient. A reclassification of non-Newtonian fluid based on the fractional constitutive equation (1) was given by (Sun et al., 2018). In the aforementioned study, the case < α < 1 was proposed to describe shear thickening fluids, 2 > α > 1 model shear thinning fluids and α = 1, viscous fluids. The partially coupled system of equations describing the transport processes is given as (Chen et al., 2015;Khan and Pop, 2010) ∂u ∂x subject to the boundary conditions v = 0, u = U0, T = Tw, C = Cw at y = 0, u = 0, T = T∞, C = C∞ as y → ∞.
Here C 0 D α y is the Caputo fractional operator of order α of the function u(y) : (0, ∞) → R, defined as The Caputo fractional operator is the most popular in applications due to the ease in imposing conditions on the fractional differential equations. Other fractional operators include the Riemann-Liouville, Riesz and the Canavati fractional derivatives. We note that the Caputo fractional operator is a regularized form of the Riemann-Liouville fractional derivative. The definitions and properties of the Caputo derivative and the aforementioned arbitrary order derivatives can be found in Miller and Ross (1993); Podlubny (1998); Atanackovic et al. (2014). We define u and v (m/s) as the horizontal and vertical velocity components respectively, να(m 1+α /s) = µα/ρ b is the fractional kinematic viscosity, ρ b and ρp (kg/m 3 ) are constant basefluid and nanoparticle densities respectively. If we assume the reference length L and velocity U0, we have the following dimensionless variables (Reα, M 2 , Le, P r, N bt, Sc) = U0L α να , where N b and N t are respectively the Brownian diffusion and thermophoresis parameters, Reα is the fractional Reynolds number, M is the magnetic parameter, Le is the Lewis number, P r is the Prandtl number, and Sc is the Schmidt number (Dhanai et al., 2015). Introducing the stream function ψ defined by U = ∂ψ/∂Y and V = −∂ψ/∂X, Eqs.
(3) to (5) become ∂ψ ∂Y We perform a scaling transformation by introducing the new variables x, y,ψ,θ,φ defined as Here, λ is the scaling parameter and ai for i = 1 . . . 5 are constants to be determined such that Eqs. (9) to (11) remain invariant under this group transformation. Upon substituting Eq. (12) into Eqs. (9) to (11), and dropping the tilde for convenience, we get λ a 3 +a 5 −a 2 −a 1 Sc ∂ψ ∂y For the system of equations to remain invariant, we set a1 = a3 = a (arbitrary) and a2 = a4 = a5 = 0. The characteristics equations associated with the transformation is given as which give the similarity functions Using the Caputo fractional operator to obtain the fractional derivative of the similarity function ψ(x, y) = xf (η) and substituting Eq. (17) into Eqs. (9) to (11) gives the following similarity equations with boundary conditions We note here that setting α = 1, Reα = 1 and M = 0, the system of equations reduces to the problem studied in Khan and Pop (2010).

NUMERICAL APPROXIMATION
Existing numerical methods for integer order differential equations have been extended to fractional order differential equations. Numerical and approximate methods such as the differential transform, homotopy, Adomian decomposition and finite difference methods have been used to solve fractional differential equations (Wang, 2006;Diethelm et al., 2002;Erjaee et al., 2011;El-Sayed et al., 2010;Momani and Odibat, 2007a,b). Despite the fact that some of these methods are considered accurate and reliable, they have certain limitations. For example, due to the history dependence associated with fractional derivatives, approximations can be computationally expensive. Fractional derivatives are defined globally, hence numerical methods that discretize non-locally would be appropriate in approximating the derivative. One such method is the spectral method. Spectral methods are synonymous with exponential rate of convergence. We discretized the dependent variables and their arbitrary order derivatives using the spectral collocation methods. The Chebyshev spectral collocation method uses trial functions to represent the functions as truncated series expansions. The trial functions are usually orthogonal basis functions such as the Lagrange, Jacobi, Chebyshev and other polynomials. A crucial difference between spectral methods and other methods is that the test functions for spectral methods are inherently continuous and infinitely differentiable global functions. Some recent studies have demonstrated the use of the operational matrix for certain orthogonal polynomials. Ahmadi Darani and Saadatmandi (2017)  usually lead to accurate solutions. However, dealing with nonlinearity and coupling can be time consuming. Hence, in this study, we follow a different approach by using approximations based on the shifted Chebyshev polynomials of the first kind and integration on the Gauss-Lobatto quadrature to obtain the coefficients of expansion. The series form of shifted Chebyshev polynomials T l,n of degree n > 0, where l is the truncation of the positive half domain, is defined (Abramowitz and Stegun, 1965;Doha et al., 2011;Ahmadi Darani and Saadatmandi, 2017) as and satisfies the orthogonality condition We discretized the dependent variables using the shifted Chebyshev Gauss Lobatto grid points η ∈ [0, l] The weight function w l (η) for the shifted Chebyshev polynomials is given by 1/ lη − η 2 and hn = cnπ/2, with c0 = 2 and cn = 1 for n ≥ 1.
If F (η) is approximated using the shifted Chebyshev polynomials and evaluated at the shifted Gauss-Lobatto collocation points, then any arbitrary order derivative of F (η) is given as Here, n,k T l,k (ηp), j, p = 0, 1, . . . , N.
and D (α) n,k = n n j=0 (−1) n−j (n + j − 1)!2 2j (n − j)!(2j)!l j Γ(j + 1) Global Digital Central ISSN: 2151-8629 We note here that if α = 1, Eq. (27) is equivalent to n cn T l,n (ηj)S l,n−1 (η k ) F (ηj), where S l,n is the shifted second kind Chebyshev polynomials of degree n ≥ 1. In order to approximate f (η) in terms of the Chebyshev polynomials, it has been assumed that f (η) is a smooth function, so it is safe to take advantage of the semi-group property of the Caputo fractional operator, therefore when approximating f α+2 (η), we can have (Matlob and Jamali, 2017;Dabiri and Butcher, 2018)

Error bound estimation
If F (η) is a square integrable function and w l (η) is a Lebesque integrable function defined in the interval (0, l), we can define a L 2 w l space in which F (η) is measurable and the norm of F (η) associated with the space is defined as and this norm is induced by the inner product If we consider TN , the space of all Chebyshev polynomials of degree ≤ N , such that TN F (η) is the approximation of F (η) on the shifted Chebyshev-Gauss-Lobatto nodes. Assume that the derivative F N +1 (η) exist and is continuous on the interval (0, l), and using the error term of the generalized Taylor's approximation of F (η), for any η, the error bound is defined as where QN = max 0<η≤l F N +1 (η) 2 . For sufficiently large N , we can see that εN w l → 0.

NUMERICAL RESULTS
In this section, we obtain approximate solutions of the equations describing heat and mass transfer in a nanofluid with space fractional constitutive relation using the technique described in Section 3. The aim is to obtain approximate solutions of the system of differential Eqs. (18)  . Some results obtained using the above scheme are compared with results in these studies. The nonlinear system is first linearized using quasilinearization as in Bellman and Kalaba (1965); Motsa et al. (2011Motsa et al. ( , 2014 [see Appendix A for the discretization of the equations]. In developing numerical schemes, a pertinent question often asked centers around the accuracy of the numerical scheme. For an approximation method, it is expected that when the approximate solutions and their derivatives are substituted in Eqs. (18) to (20), the residual errors should vanish. That is, given the points ηj ∈ [0, l] for j = 0, . . . , N , we should expect and so the skin friction is given as      Table 5 Comparison of the heat and mass transfer coefficients at the wall for different values of N t and N b with α = 1 and fixed values for P r = Le = 10, Reα = 1 and M = 0 Present Study Khan and Pop (2010) Present Study Khan and Pop (2010) Table 6 Numerical values of the shear stress at the wall, rate of heat and mass transfer coefficients for different values of the fractional order α, magnetic parameter M , and Prandtl number P r when Le = 5, N t = N b = 0.5. To establish the accuracy of the numerical scheme, we compare Equation (42) with the approximate solutions. Table 4 shows the comparison between the exact and approximate values of the skin friction coef-  Table 4 are obtained with N = 10 and l = 5.9. The heat and mass transfer coefficients {θ (0), φ (0)} obtained using this method are compared with earlier results by Khan and Pop (2010). Table 5 shows that the new results are in agreement with the published data. Figure 2 shows that the solutions converge for different number of collocation points, especially in the equation where we have the fractional order. In Table 6, the approximate values of the skin friction and heat and mass transfer rates are presented. Figs. 3 and 4 respectively show the variation in the heat and mass transfer rates with the Brownian diffusion and thermophoresis parameters for different values of the fractional order α. It can be clearly seen that the heat and mass transfer rates increase as α increases. Fig. 3(a) shows that the heat transfer rate decreases as the mass diffusion parameter N b increases, while Fig. 3(b) shows that the mass transfer rate is a monotonically increasing function of N b. From Fig. 4, we can see that the mass transfer rate decreases as the thermal diffusion parameter increases, while the heat transfer rate strictly decreases. Previous studies have suggested that increasing the number of collocation points for spectral methods does not necessarily improves the accuracy (Dabiri and Butcher, 2017b,a; Zayernouri and Karniadakis, 2014).

CONCLUSION
This study proposed a new numerical technique for finding solutions of fractional differential equations. The method uses Chebyshev polynomials as basis function for the series solution and interpolating using Gauss-Lobatto quadrature. Differentiation matrices for fractional and integer order derivatives were introduced. To test its accuracy, the method was used to solve the model equations that describe momentum, heat and mass transfer in a nanofluid with fractional constitutive relation. The residual error norms in the solutions were found to be less or equal to 10 −13 , thus demonstrating the accuracy and convergence of the method.