OPEN ACCESS
We focused on developing an accurate numerical scheme for the flow of a fractionalorder Oldroyd–B fluid model with the nonisothermal property. In many cases, the direct application of the Chebyshev tau method using the operational matrix of the Chebyshev polynomials usually leads to an accurate solution. However, in some cases, dealing with non–linearity and coupling can be tedious. In this study, we present a numerical method based on Chebyshev polynomials of the first kind and interpolation using Gauss–Lobatto quadrature. The coefficients of the series expansion of the pseudospectral method are obtained through integration of the Chebyshev polynomials orthogonality condition. The numerical results show that the scheme is accurate and reliable. The effects of the fractionalorder objective stress rate of the Oldroyd–B fluid on the velocity and shear stress are also presented. The error bound theorems presented in this study support the findings of the numerical computations.
MHD fluid, nonisothermal flow, fractional calculus, Chebyshev – Gauss – Lobatto quadrature, fractional Oldroyd–B fluid
The Oldroyd–B fluid is a rate type model for viscoelastic fluids that have been studied extensively since its formulation by Oldroyd [1]. This model is part of a large class of rate type viscoelastic fluids models, which include, but are not limited to the upper convected Maxwell model, Jeffrey model and Burgers’ model. The OldroydB model gives a good representation of the rheological response of viscoelastic fluids in shear flow. The Cauchy stress tensor for the classical Oldroyd–B constitutive model has the form
$S_{i j}=p \delta_{i j}+\tau_{i j}$ (1)
$[1+\lambda \overset{\nabla }D] \tau_{i j}=\mu\left[1+\lambda_{r} \overset{\nabla }{D}\right]e_{ij}$ (2)
where, τ_{ij} is the extra stress tensor, e_{ij} is the stress tensor. The material parameters μ, λ and λ_{r} are the viscosity, stress relaxation time and retardation time respectively and $\overset{\nabla }{D}$ is the upper convected time derivative. The constitutive model of the Oldroyd–B is a generalization of the Maxwell fluid model [2]. Unlike the Maxwell model, the Oldroyd–B fluid uses the objective stress rate that takes into account frame indifference of the deformation rate. The deformation rate is represented by the upper convected derivative in Eq. (2). Because of the memory retention characteristics of viscoelastic fluids, recent studies have suggested using non–local time derivatives [3]. Non–local objective stress rate exhibits complex dynamical and viscoelastic behaviours. Taking advantage of the inherent non–locality and memory retention characteristics of fractional derivatives, we generalize the Oldroyd–B constitutive relation by replacing the local time derivatives in Eq. (2) with fractional time derivative [4]:
$\left[1+\lambda^{\alpha}\overset{\nabla }{{\mathrm{D}}^{\alpha}}\right] \tau_{i j}=\mu\left[1+\lambda_{r}^{\beta}\overset{\nabla }{{\mathrm{D}}^{\beta}}\right] e_{i j}$ (3)
where,
$\overset{\nabla }{D^{\alpha}}=\frac{\partial^{\alpha} \tau}{\partial t^{\alpha}}+\mathbf{U} \cdot \nabla \tau\tau \nabla \mathbf{U}(\nabla \mathbf{U})^{\mathrm{T}} \tau$ (4)
In the constitutive model Eq. (3), 0<α, β≤1 and μ, λ, λ_{r}>0. If α>β, it should be noted that the model is physically unrealistic. This case corresponds to an increasing relaxation function [5]. Hence, the constraint on the order of the derivatives must satisfy 0<α≤β≤1. It is evident that if α=β=1, the constitutive relation is the classical Oldroyd–B fluid. Setting λ_{r} to zero in Eq. (3) corresponds to the fractional Maxwell fluid, if λ=0, the model is equivalent to the fractional secondgrade fluid, if 0<λ_{r}<λ, we obtain the fractional Jeffrey fluid [6] and when λ=λ_{r}=0, the model corresponds to the classical Newtonian fluid. If we consider an incompressible fluid with constant pressure and using established notation, the momentum conservation equation for an unsteady flow of a magnetohydrodynamic (MHD) fluid is:
$\rho\left[\frac{\partial \mathbf{U}}{\partial t}+\mathbf{U} \cdot \nabla \mathbf{U}\right]=\nabla \tau_{i j}\sigma B_{0}^{2} \mathbf{U}$ (5)
where, ρ is the fluid density, σ is the electrical conductivity of the fluid and B_{0} is the applied magnetic field. To write the conservation equation and the shear stress relation in non–dimensional form, we assume a reference length L and velocity $\mathbf{U}_{0},$ so that we define $t=\tilde{t} \mathbf{L} \mathbf{U}_{0}^{1}, \mathbf{y}=\tilde{\mathbf{y}} \mathbf{L}$ for all $\mathbf{y}$ in the domain of the flow, $\mathbf{U}(\mathbf{y}, t)=\mathbf{U}_{0} \widetilde{\mathbf{U}}(\tilde{\mathbf{y}}, \tilde{t}), \tau_{i j}(\mathbf{y}, t)=$$\mathbf{U}_{0} \mu \mathbf{L}^{1} \tilde{\tau}_{i j}(\tilde{\mathbf{y}}, \tilde{t}),$ and $\lambda_{r}=\tilde{\lambda}_{r} \mathbf{L} \mathbf{U}_{0}^{1} .$ Therefore, we write Eq. (3) and Eq. (5) in the dimensionless form (dropping the tilde for convenience):
$\operatorname{Re}\left[\frac{\partial \mathbf{U}}{\partial t}+\mathbf{U} \cdot \nabla \mathbf{U}\right]=\nabla \tau_{i j}H a^{2} \mathbf{U}$ (6)
$\left[1+W e^{\alpha}\overset{\nabla }{{D}^{\alpha}}\right] \tau_{i j}=\left[1+\lambda_{r}^{\beta}\overset{\nabla }{{\mathrm{D}}^{\beta}}\right] e_{i j}$ (7)
Here, $R e=\rho \mathbf{L U}_{0} \mu^{1}$ is the Reynolds number, $H a=$ $B_{0} \mathbf{L}(\sigma / \mu)^{1 / 2}$ is the Hartmann number, $W e=\mathbf{U}_{0} \lambda \mathbf{L}^{1}$ is the Weissenberg number, λ_{r} is a dimensionless retardation parameter. Some special cases of interest include the Jeffrey fluid (0<λ_{r}<We), Maxwell fluid (λ_{r}=0), secondorder fluid (We=0), Newtonian fluid (We=λ_{r}=0) and hydrodynamic fluids (Ha=0).
Exact solutions to fractional partial differential equations are, in general, challenging to find. Even when found, the solution may contain complicated integrals and special functions that have to be approximated numerically. Several studies have proposed various methods, both numerical and analytic, for solving fractional differential equations [712]. However, recent studies have identified spectral methods as efficient for approximating solutions of fractional partial differential equations. Spectral methods are favoured because of their spectral rate of convergence and accuracy, especially for differential equations with sufficiently smooth solutions. Doha et al. [13] and Atabakzadeh et al. [14] presented the operational matrix of the shifted Chebyshev polynomial and used the polynomial to approximate multiorder fractional ordinary differential equations. Liu et al. [15] presented numerical solutions to multiterm variableorder fractional ordinary differential equations using the Chebyshev polynomial of the second kind as the basis function. Operational matrices of the shifted Chebyshev wavelets, generalized Laguerre polynomials, and Bernstein polynomials are presented in the studies of Benattia and Belghaba [16], Bhrawy and Alghamdi [17], Baseri et al. [18] respectively. These matrices are used to approximate the solutions of fractional differential equations. The numerical results that were reported in the studies as mentioned earlier are typical of spectral method, and they exhibit the spectral convergence property of spectral methods. From literature, many authors have used the spectral tau method to solve fractional differential equations with different orthogonal polynomials as basis functions. In this study, we propose a spectral method that uses the shifted Chebyshev polynomials of the first kind as basis functions and interpolates using Gauss–Lobatto quadrature. We obtain the coefficients of the series expansion in the spectral method by integrating the orthogonal condition of the Chebyshev polynomials using ChebyshevGaussLobatto quadrature. We then apply the resulting fractional differentiation matrix to solve fractional partial differential equations that describe the unsteady boundary layer flow of a generalized MHD Oldroyd–B fluid and a non–isothermal flow of an Oldroyd–B fluid (see Section 4). We also investigate the effects of derivatives of arbitrary order in Eq. (6) and Eq. (7) on fluid velocity and shear stress.
The most used definitions of fractional operators are the RiemannLiouville and Caputo fractional operators. In the RiemannLiouville case, we have the definition
${ }^{R L} \mathrm{D}_{y}^{\alpha} H(y)=\frac{1}{\Gamma(n\alpha)} \frac{d^{n}}{d y^{n}} \int_{0}^{y} \frac{H(\varsigma)}{(y\varsigma)^{\alpha+1n}} d \varsigma$ (8)
and the Caputo case is defined as
$\begin{array}{ll} & ^C_0 D_{y}^{\alpha} H(y) \\ & =\left\{\begin{array}{cl}\frac{1}{\Gamma(n\alpha)} \int_{0}^{y} \frac{H^{(n)}(\zeta)}{(y\zeta)^{\alpha+1n}} d \zeta, & & n1<\alpha<n \\ \frac{d^{n} H}{d y^{n}}, & \alpha=n\end{array}\right.\end{array}$ (9)
In both cases, n1<α<n, and $n \in \mathbb{N}$ and H(y) is a continuously bounded function with n derivatives in [0, L], for χ>0. For an indepth exploration of these definitions, see Podlubny [19], Atangana [20]. Using the Caputo operator, the following holds
${ }^{C} D^{\alpha} K=0, \quad \mathrm{K}$ is a constant (10)
$c_{D^{\alpha} y^{j}}=\left\{\begin{array}{cc}0 & \text { for } j \in \mathbb{N}_{0} \text { and } j<\lceil\alpha\rceil \\ \frac{\Gamma(j+1)}{\Gamma(j+1\alpha)} y^{j\alpha} & \text { for } j \in \mathbb{N}_{0}, j \geq[\alpha]\end{array}\right.$ (11)
Considering that most physical problems are defined in [0, L], L being a truncation of the semiinfinite domain, we define the series form of the shifted Chebyshev polynomials T_{L,n} of degree n>0 as ([21, 22])
$T_{L, n}=n \sum_{j=0}^{n} \frac{(1)^{nj}(n+j1) ! 2^{2 j}}{(nj) !(2 j) ! L^{j}} y^{j}$ (12)
which satisfies the orthogonality condition
$\int_{0}^{L} T_{L, n}(y) T_{L, m}(y) w_{L}(y) d y=\delta_{m n} h_{n}$ (13)
Here, the weight function of the shifted Chebyshev polynomials is given by $w_{L}(y)=1 / \sqrt{L yy^{2}}$ and h_{n}=c_{n} π/2, with c_{0}=2 and c_{n}=1 for n≥1.
At this point, we discuss some lemmas and theorems that are used in developing the approximation of the fractionalorder derivative of a squareintegrable function that is expanded by a shifted Chebyshev series and integrated using the ChebyshevGaussLobatto quadrature.
Remark 2.1: For the shifted Chebyshev–Gauss–Lobatto quadrature, the Christoffel numbers are the same as those of the ChebyshevGaussLobatto quadrature. The shifted Gauss–Lobatto nodes are defined as ([21])
$y_{j}=\frac{L}{2} \cos \left(\frac{\pi j}{N}\right)+\frac{L}{2}$ (14)
and the associated Christoffel weight number w_{L,j}=π/c_{i }N, 0≤j≤N, where c_{0}=c_{N}=2 and c_{j}=1 for j=1,2,…, N1.
Lemma 2.2: Assume that u(y) is an integrable function defined on the domain[0,L], then the function can be expanded in terms of shifted Chebyshev polynomials as a N+1 truncated series:
$u_{N}(y)=\sum_{n=0}^{N} u_{n} T_{L, n}(y)$ (15)
where, the coefficients u_{n} satisfy the orthogonality condition, which is given in discrete form as
$u_{n}=\frac{1}{h_{n}} \sum_{j=0}^{N} \frac{\pi}{c_{j} N} u\left(y_{j}\right) T_{L, n}\left(y_{j}\right), \quad n=0, \ldots, N$ (16)
Lemma 2.3: Let T_{L,n} (y) be a n–th order shifted Chebyshev polynomials, then the a–th order derivative based on the Caputo’s definition is defined as [13, 14]
$\mathrm{D}^{\alpha} T_{L, n}(y)=\sum_{k=0}^{N} \mathrm{D}_{n, k}^{(\alpha)} T_{L, k}(y)$ (17)
where,
$D_{n, k}^{(\alpha)}$$=n \sum_{j=0}^{n} \frac{(1)^{nj}(n+j1) ! 2^{2 j}}{(nj) !(2 j) ! L^{j}} \frac{\Gamma(j+1)}{\Gamma(j\alpha+1)} q_{j, k}$ (18)
and q_{j}_{,k} are entries of a matrix and these entries are defined as (Ahmadi  Darani and Saadatmandi [22], Ahmadi  Darani and Nasiri [23])
$q_{j, k}=\left\{\begin{array}{ll}0 & j=0,1, \ldots,[\alpha]1 \\ \frac{k \sqrt{\pi}}{h_{k}} \sum_{r=0}^{k} \frac{(1)^{kr}(k+r1) ! 2^{2 r}}{(kr) !(2 r) !} L^{j\alpha} \frac{\Gamma\left(j\alpha+r+\frac{1}{2}\right)}{\Gamma(j\alpha+r+1)}, & j=\lceil\alpha\rceil,\lceil\alpha\rceil+1, \ldots, N \\ & k=0,1 \ldots N\end{array}\right.$ (19)
The proof of this lemma can be found in Atabakzadeh et al. [14] and Doha et al. [13].
Theorem 2.4: Assume u(y) is a continuously bounded and integrable function defined on the truncated semiinfinite domain [0,L]. If u(y) is approximated using the shifted Chebyshev polynomials and evaluated at the shifted GaussLobatto collocation points, then any arbitrary order derivative of u(y) is given by
$\mathrm{D}^{\alpha} u_{N}(y)=\sum_{j=0}^{N} \mathrm{D}_{j, p}^{\alpha} u\left(y_{j}\right)$ (20)
Here,
$\begin{aligned} D_{j, p}^{\alpha}=\frac{\pi}{c_{j} N} & \sum_{n=0}^{N} \sum_{k=0}^{N} \frac{1}{h_{n}} T_{L, n}\left(y_{j}\right) \mathrm{D}_{n, k}^{(\alpha)} T_{L, k}\left(y_{p}\right) \\ & j p=0,1, \ldots, N \end{aligned}$ (21)
Proof. If we consider the first N+1 shifted Chebyshev polynomials, and a combination of the results of Lemma 2.2 and Lemma 2.3, the proof is completed.
Hitherto, we have only focused on functions of one variable. We now focus on functions of two variables, precisely time and spatial variables. We will use the results of the lemmas and theorem stated above to obtain the approximated values of a bivariate function and its fractionalorder partial derivatives. The most straightforward approach to expanding a function of two variables is to take the tensor product of the expansion in each variable [24]. For a function u(y,t), the expansion by shifted Chebyshev polynomials is given as
$u_{N_{y}, N_{t}}(y, t)=\sum_{n=0}^{N_{y}} \sum_{m=0}^{N_{t}} u_{n, m} T_{L, n}(y) T_{\mathrm{T}, m}(t)$ (22)
where, u_{n,m} are entries of the matrix which consist of the coefficients of expansion.
Proposition 2.5: If u(y,t) is a function defined on the rectangular domain [0,L]×[0,T], then it can be approximated at the shifted GaussLobatto nodes in terms of N+1 shifted Chebyshev polynomials so that we have
$\operatorname{vec}(\vec{u})=\left[I_{t} \otimes I_{y}\right] \operatorname{vec}(\vec{u})$ (23)
Here, both I_{y} and I_{t} are interpolation matrices (Lemma 2.2) in y and t respectively, and $\operatorname{vec}(\vec{u})$ is formed by piling columns of {u(y_{i},t_{j})  i=0,…, N_{y}, j=0,…,N_{t}} into a vector.
Proposition 2.6: Any arbitrary order, say α>0 derivatives of u(y,t) can also be expanded and approximated in terms of the shifted Chebyshev polynomials and shifted Gauss–Lobatto nodes. The α–th order derivatives in terms of both variables are given as
$\operatorname{vec}\left(\vec{u}_{\alpha}\right)=\left[I_{t} \otimes \mathrm{D}_{y}^{\alpha}\right] \operatorname{vec}(\vec{u})$
$\operatorname{vec}\left(\vec{u}^{\alpha}\right)=\left[\mathrm{D}_{t}^{\alpha} \otimes I_{y}\right] \operatorname{vec}(\vec{u})$
$\operatorname{vec}\left(\vec{u}_{\alpha}^{\alpha}\right)=\left[\mathrm{D}_{t}^{\alpha} \otimes \mathrm{D}_{y}^{\alpha}\right] \operatorname{vec}(\vec{u})$
where, the superscript indicates the derivative with respect to the temporal variable and the subscript is the derivative with respect to the spatial variable and D^{α} is from Theorem 2.4.
Assume that $u(y)$ is a squareintegrable function and $w_{L}(y)$ is a Lebesgue integrable function defined in the interval $\mathbb{I}=$ $[0, L] .$ Then, we can define a $\mathbf{L}_{w_{L}}^{2}$ space in which $u(y)$ is measurable and the norm $\u(y)\_{w_{I} .}$ is defined as
$\u(y)\_{w_{L}}=\left(\int_{\mathbb{I}}u(y)^{2} w_{L}(y) d y\right)^{1 / 2}<\infty$ (24)
such that the norm is induced by the inner product
$\langle u(y), \hat{u}(y)\rangle=\int_{0}^{L} u(y) \hat{u}(y) w_{L}(y) d y$ (25)
If u(y) is the exact solution and u_{N} (y) is the approximated solution based on the Chebyshev polynomials given in Eq. (15) with the coefficients u_{n} obtained by the normalized inner product
$u_{n}=\frac{\left\langle u_{N}(y), T_{L, n}(y)\right\rangle}{\left\T_{L, n}(y)\right\_{w_{L}}^{2}}$ (26)
We can define an error bound for the approximation in the $\mathbf{L}_{w_{L}}^{2}$ norm.
Theorem 3.1 (Error estimation for a single variable approximation): Given the shifted interpolation nodes defined in Eq. (14) and let P_{N} u(y) be the approximation through these nodes given in Eq. (15), where P_{N }is the space of all Chebyshev polynomials of degree less than or equal to N. Assume that d^{N}^{+1}u/dy^{N}^{+1} exist and is continuous on the interval $\mathbb{I}$, then the error bound is defined as
$\left\u(y)\mathbf{P}_{N} u(y)\right\$$\leq \frac{1}{(\Gamma(N+2))^{2}}\left(\max _{0<y \leq L}\left\frac{d^{N+1} u(y)}{d y^{N+1}}\right\right)^{2} L^{2 N+2} \sqrt{\pi} \frac{\Gamma\left(2 N+\frac{5}{2}\right)}{\Gamma(2 N+3)}$ (27)
Proof. Consider the generalized Taylor’s approximation of u(y) in which the error bound is known as
$\leftR_{N}(y)\right \leq \frac{y^{N+1}}{\Gamma(N+2)} \max _{0<y \leq L}\left\frac{d^{N+1} u(y)}{d y^{N+1}}\right$ (28)
then for any y in the collocation points
$\left\u(y)\mathbf{P}_{N} u(y)\right\_{w_{L}}^{2} \leq\left\R_{N}(y)\right\_{w_{L}}^{2}$$\leq \frac{1}{(\Gamma(N+2))^{2}}\left(\max _{0<y \leq L}\left\frac{d^{N+1} u(y)}{d y^{N+1}}\right\right)^{2} \int_{0}^{L} \frac{y^{2 N+2}}{\sqrt{L yy^{2}}} d y$ (29)
$=\frac{1}{(\Gamma(N+2))^{2}}\left(\max _{0<y \leq L}\left\frac{d^{N+1} u(y)}{d y^{N+1}}\right\right)^{2} L^{2 N+2} \sqrt{\pi} \frac{\Gamma\left(2 N+\frac{5}{2}\right)}{\Gamma(2 N+3)}$ (30)
This leads to the desired result.
Theorem 3.2: Let $u: \mathbb{I} \times \mathbb{J} \rightarrow \mathbb{R}$ be a continuously differentiable function such that at least (N_{y}+1)th partial derivative with respect to y, (N_{t}+1)th partial derivative with respect to t, and (N_{y}+N_{t}+2)th mixed derivative with respect to y and t exist, then based on the mean value theorem, the following remainder formula holds ([25])
$\begin{aligned}\leftR_{N_{y}, N_{t}}(y, t)\right \leq & \frac{K_{1} y^{N_{y}+1}}{\Gamma\left(N_{y}+2\right)}+\frac{K_{2} t^{N_{t}+1}}{\Gamma\left(N_{t}+2\right)} \\ &+\frac{K_{3} y^{N_{y}+1} t^{N_{t}+1}}{\Gamma\left(N_{y}+2\right) \Gamma\left(N_{t}+2\right)^{\prime}} \end{aligned}$ (31)
where, $\mathbb{I}=\lceil 0, L\rceil, \mathbb{J}=\lceil 0, \mathrm{T}\rceil$ and K_{1}, K_{2}, K_{3} are constants, defined as
$K_{1}=\sup \left\{\left\frac{\partial^{{N} y^{+1}} u(y, t)}{\partial y^{N_{y}+1}}\right: y, t \in \mathbb{I} \times \mathbb{J}\right\}$
$K_{2}=\sup \left\{\left\frac{\partial^{N_{t}+1} u(y, t)}{\partial t^{N_{t}+1}}\right: y, t \in \mathbb{I} \times \mathbb{J}\right\}$ (32)
$K_{3}=\sup \left\{\left\frac{\partial^{N_{y}+N_{t}+2} u(y, t)}{\partial y^{N_{y}+1} \partial t^{N_{t}+1}}\right: y, t \in \mathbb{I} \times \mathbb{I}\right\}$ (33)
respectively.
Theorem 3.3 (Error bound for functions of 1+1 variables approximation): We define the error bound for the approximation of a function of two variables as
$\left\R_{N_{y}, N_{t}}(y, t)\right\_{w_{L}, w_{\mathrm{T}}}$$\leq \frac{K_{1}^{2} \pi \sqrt{\pi} L^{2 N_{y}+2} \Gamma\left(2 N_{y}+\frac{5}{2}\right)}{\left(\Gamma\left(N_{y}+2\right)\right)^{2}\Gamma\left(2 N_{y}+3\right)}+ \frac{K_{2}^{2} \pi \sqrt{\pi} \mathrm{T}^{2 N_{t}+2} \Gamma\left(2 N_{t}+\frac{5}{2}\right)}{\left(\Gamma\left(N_{t}+2\right)\right)^{2}\Gamma\left(2 N_{t}+3\right)}$
$+\frac{K_{3}^{2} \pi L^{2 N_{y}+2} T^{2 N_{t}+2}}{\left(\Gamma\left(N_{y}+2\right)\right)^{2}\left(\Gamma\left(N_{t}+2\right)\right)^{2}} \frac{\Gamma\left(2 N_{t}+\frac{5}{2}\right) }{\Gamma\left(2 N_{t}+3\right)}\frac{\Gamma(2N_y +\frac{5}{2})}{\Gamma\left(2 N_{y}+3\right)}$ (34)
Proof. Given that $\mathbf{P}_{N_{y}, N_{t}} u(y, t)$ is the space of all Chebyshev polynomials which approximate u(y,t). In a similar sense as in Theorem 3.1, for any points y,t in the collocation points and using the relation in Eq. (31), we define
$\left\u(y, t)\mathbf{P}_{N_{y}, N_{t}} u(y, t)\right\_{w_{L}, w_{\mathrm{T}}}^{2}$$\leq\left\R_{N_{y}, N_{t}}(y, t)\right\_{w_{L}, w_{\mathrm{T}}}^{2}$ (35)
$\leq \int_{\mathbb{I} \times \mathbb{J}}\left(\left\frac{K_{1} y^{N_{y}+1}}{\Gamma\left(N_{y}+2\right)}\right^{2}+\left\frac{K_{2} t^{N_{t}+1}}{\Gamma\left(N_{t}+2\right)}\right^{2}\right.$
$\left.+\left\frac{K_{3} y^{N_{y}+1} t^{N_{t}+1}}{\Gamma\left(N_{y}+2\right) \Gamma\left(N_{t}+2\right)}\right^{2}\right) \frac{1}{\sqrt{\mathrm{T} tt^{2}}} \frac{1}{\sqrt{L yy^{2}}} d y d t$ (36)
$=\frac{K_{1}^{2}}{\left(\Gamma\left(N_{y}+2\right)\right)^{2}} \int_{\mathbb{I} \times \mathbb{J}} \frac{1}{\sqrt{\mathrm{T} tt^{2}}} \frac{y^{2 N_{y}+2}}{\sqrt{L yy^{2}}} d y d t$
$+\frac{K_{2}^{2}}{\left(\Gamma\left(N_{t}+2\right)\right)^{2}} \int_{\mathbb{I} \times \mathbb{J}} \frac{t^{2 N_{t}+2}}{\sqrt{\mathrm{T} tt^{2}}} \frac{1}{\sqrt{L yy^{2}}} d y d t$
$+\frac{K_{3}^{2}}{\left(\Gamma\left(N_{t}+2\right)\right)^{2}\left(\Gamma\left(N_{y}+2\right)\right)^{2}} \int_{\mathbb{I} \times \mathbb{J}} \frac{t^{2 N_{t}+2}}{\sqrt{T tt^{2}}} \frac{y^{2 N_{y}+2}}{\sqrt{L yy^{2}}} d y d t$ (37)
The integrals in Eq. (37) are evaluated as
$\int_{\mathbb{I} \times \mathbb{J}} \frac{1}{\sqrt{\mathrm{T} tt^{2}}} \frac{y^{2 N_{y}+2}}{\sqrt{L yy^{2}}} d y d t$
$=\pi \sqrt{\pi} L^{2 N_{y}+2} \frac{\Gamma\left(2 N_{y}+\frac{5}{2}\right)}{\Gamma\left(2 N_{y}+3\right)}$ (38)
$\int_{\mathbb{I} \times \mathbb{J}} \frac{t^{2 N_{t}+2}}{\sqrt{\mathrm{T} tt^{2}}} \frac{1}{\sqrt{L yy^{2}}} d y d t$
$=\pi \sqrt{\pi} \mathrm{T}^{2 N_{t}+2} \frac{\Gamma\left(2 N_{t}+\frac{5}{2}\right)}{\Gamma\left(2 N_{t}+3\right)}$ (39)
$\int_{\mathbb{I} \times \mathbb{J}} \frac{t^{2 N_{t}+2}}{\sqrt{\mathrm{T} tt^{2}}} \frac{y^{2 N_{y}+2}}{\sqrt{L yy^{2}}} d y d t$
$=\pi L^{2 N_{y}+2} T^{2 N_{t}+2} \frac{\Gamma\left(2 N_{t}+\frac{5}{2}\right)}{\Gamma\left(2 N_{t}+3\right)} \frac{\Gamma\left(2 N_{y}+\frac{5}{2}\right)}{\Gamma\left(2 N_{y}+3\right)}$ (40)
Substituting Eq. (38) to Eq. (40) into Eq. (37) completes the proof. If we consider $u_{N_{y}, N_{t}}(y, t),$ the $\mathbf{L}_{w_{L}, w_{\mathrm{T}}}^{2}$ orthogonal projection of $u(y, t)$ onto $\mathbf{P}_{N_{y}, N_{t}},$ then
$\begin{aligned} \mid u(y, t)u_{N_{y}, N_{t}}(&y, t) \mid \\ &=\mid u(y, t)\mathbf{P}_{N_{y}, N_{t}} u(y, t) \\ &+\mathbf{P}_{N_{y}, N_{t}} u(y, t)u_{N_{y}, N_{t}}(y, t) \mid \end{aligned}$ (41)
$\begin{aligned} \leq\leftu(y, t)\mathbf{P}_{N_{y}, N_{t}} u(y, t)\right & \\+\leftu_{N_{y}, N_{t}}(y, t)\mathbf{P}_{N_{y, N_{t}}} u(y, t)\right \end{aligned}$ (42)
Based on Theorem 3.3, Eq. (42) tends to zero as $\infty$.
In this section, we apply the numerical scheme using the Chebyshev–Gauss–Lobatto quadrature to the flow of an MHD and non–isothermal generalized Oldroyd–B fluid.
4.1 Flow of an MHD Oldroyd–B fluid
We investigate an unsteady unidirectional flow of an MHD OldroydB fluid occupying the upper halfplane bounded by an impermeable wall at y=0. The y–axis is perpendicular to the rigid flat plate, and the plate oscillates parallel to itself. The fluid is assumed to be at rest until startup. The momentum and shear stress equations for this unsteady flow are given as [26]
$\operatorname{Re}\left[\frac{\partial}{\partial t}+W e^{\alpha} \frac{\partial^{\alpha+1}}{\partial t^{\alpha+1}}\right] u\left[\frac{\partial^{2}}{\partial y^{2}}+\lambda_{r}^{\beta} \frac{\partial^{\beta+2}}{\partial t^{\beta} \partial y^{2}}\right] u$
$+H a^{2}\left[1+W e^{\alpha} \frac{\partial^{\alpha}}{\partial t^{\alpha}}\right] u=0$ (43)
$\left[1+W e^{\alpha} \frac{\partial^{\alpha}}{\partial t^{\alpha}}\right] \tau=\left[\frac{\partial}{\partial y}+\lambda_{r}^{\beta} \frac{\partial^{\beta+1}}{\partial t^{\beta} \partial y}\right] u$ (44)
with initial conditions
$u(y, 0)=\frac{\partial}{\partial t} u(y, 0)=0, \quad \tau(y, 0)=0$ (45)
and boundary conditions with a sine oscillation at the wall and dimensionless frequency of oscillation ω
$u(0, t)=\sin (\omega t), \quad u(y, t)=0$ as $y \rightarrow \infty, \quad t>0$ (46)
We seek a solution in the form of truncated Chebyshev polynomials
$\begin{aligned} u(y, t) &=\sum_{n=0}^{N_{y}} \sum_{m=0 }^{N_{t}} u_{n, m} T_{L, n}(y) T_{\mathrm{T}, m}(t) \\ \tau(y, t) &=\sum_{n=0}^{N_{y}} \sum_{m=0}^{N_t} \tau_{n, m} T_{L, n}(y) T_{\mathrm{T}, m}(t) \end{aligned}$ (47)
where, the coefficients $\hat{u}_{n, m}$ and $\hat{\tau}_{n, m}$ are defined by Lemma 2.2. Using the results of Proposition 2.6, we resolve Eq. (43) and Eq. (44) into a system of discrete equations
$\begin{aligned} \operatorname{Re}\left[\left[\left(\mathrm{D}_{t}^{1} \otimes I_{y}\right)+W e^{\alpha}\left(D_{t}^{\alpha+1} \otimes I_{y}\right)\right]\right.& \\\left[\left(I_{t} \otimes D_{y}^{2}\right)+\lambda_{r}^{\beta}\left(D_{t}^{\beta}\right.\right.& \\\left.\left.\otimes D_{y}^{2}\right)\right]+H a^{2}\left[\left(I_{t} \otimes I_{x}\right)\right.& \\\left.+W e^{\alpha}\left(D_{t}^{\alpha} I_{y}\right)\right] & v e c(\vec{u})=v e c(\overrightarrow{0}) \end{aligned}$ (48)
$\left[\left(I_{t} \otimes I_{y}\right)+W e^{\alpha}\left(D_{t}^{\alpha} \otimes I_{y}\right)\right] \operatorname{vec}(\vec{\tau})$$=\left[\left(I_{t} \otimes D_{y}^{1}\right)+\lambda_{r}^{\beta}\left(\mathrm{D}_{t}^{\beta}\right.\right.$$\left.\otimes \mathrm{D}_{y}^{1}\right) \mid \operatorname{vec}(\vec{u})$ (49)
where, $\operatorname{vec}(\overrightarrow{0})$ is a vector of zeroes of size (N_{y}+1) times (N_{t}+1) and each system is solved independently by matrix inversion.
The accuracy of the solutions for different values of α, β, N_{t}, N_{y} are presented in Table 1 in terms of the maximum residual error and condition number. The condition number was computed using the python package ‘numpy.linalg.cond’. The results in Figure 1 are shown for two frequencies of oscillation and α=0.25, β=0.75, t=10, Re=0.1 and Ha=0.5. We consider the cases: We=1<λr=1.5, Jeffrey fluid (λr=1.2<We=1.5), Maxwell fluid (We=1.5, λr=0), secondgrade fluid (We=0, λr=1.2) and Newtonian fluid. The figure shows the velocity and tangential shear stress profiles. In Figure 1a, the velocity profiles are plotted for special cases of the fractional Oldroyd–B fluid. In all cases, the oscillation of the fluid velocity decays as it approaches the boundary layer region. The fractional secondgrade fluid has the lowest amplitude of oscillation, while the fractional Maxwell fluid has the highest amplitude. This result is in agreement with the study of Jamil et al. [26]. In Figure 1b, the shear stress profiles for fixed time t are plotted for special cases. This figure shows that the case We=0 (secondgrade fluid) has the highest shear stress at the wall, while the case when λ_{r}=0 (Maxwell fluid) has the smallest shear stress. Figure 2 and Figure 3 illustrate the behaviour of the velocity and shear stress for different values of the fractional orders α and β. Figure 2 shows that as the order α increases, the shear stress decreases, while the net effect of increasing β is that the shear stress increases. In Figure 3, it can be seen that as β increases, the amplitude of oscillation becomes smaller while the amplitude becomes more significant as α increases. This result, in effect, shows that the likelihood of having a region of reverse flow close to the wall is less likely for fixed β, decreasing α and fixed α, increasing β. We remark that the degree of α is related to the relaxation time. Hence, the lack of reverse flow or small amplitude of oscillation for small values of α can be associated with the short memory of the fluid and slow response to shear stress.
Table 1. Residual and condition number for the problem in Section 4.1
(α,β) 
(N_{t},N_{y}) 
Res_{u} 
CN_{u} 
Res_{τ} 
CN_{τ} 
(0.25,0.5) 
(5,5) 
4.1482e16 
203.5 
1.3210e15 
17.35 

(10,10) 
2.8511e14 
6.0310e+04 
7.2893e15 
18.99 

(15,15) 
1.3015e07 
6.1347e+11 
1.7146e07 
8.6119e+04 
(0.5,0.75) 
(5,5 ) 
9.0899e16 
313.1 
1.1055e15 
24.02 

(10,10 ) 
3.1403e14 
1.5079e+03 
5.2698e15 
44.92 

(15,15 ) 
3.5876e07 
5.6432e+11 
1.2344e06 
4.1606e+06 
(0.75,1) 
( 5,5 ) 
5.7333e16 
477.8 
2.0091e15 
36.39 

(10,10) 
4.8934e14 
4.3557e+04 
7.8238e15 
104.4 

(15,15) 
1.2859e07 
2.2473e+11 
3.4619e07 
2.7625e+06 
(1,1) 
(5,5 ) 
7.8409e16 
625.5 
1.4658e15 
60.11 

(10,10) 
7.6247e14 
7.0962e+04 
1.7502e14 
236.8 

(15,15 ) 
1.1289e08 
2.6527e+10 
3.9879e07 
4.1076e+06 
(a)
(b)
Figure 1. Velocity and shear stress of special cases of the fractional OldroydB fluid with ω=70+nπ/2
(a)
(b)
Figure 2. Variation of the tangential shear stress profiles with ω=10+π/2 for (a) different α with β=1 and (b) different β with α=0.25
(a)
(b)
Figure 3. Velocity profiles with ω=10+π/2 for: (a) different α with β=1 and (b) different β with α=0.25
4.2 The nonisothermal flow of an Oldroyd–B fluid
Here we consider the unsteady non–isothermal flow of a generalized Oldroyd–B fluid. The thermal property of the fluid (γ), the specific heat (c) and the thermal conductivity (k) are constant and isotropic, except the viscosity which is assumed to be temperaturedependent with reference viscosity, μ_{0}. We use [27]
$\mu(T)=\frac{\mu_{0}}{1+\gamma\left(TT_{0}\right)}$ (50)
to model the nonisothermal viscosity. The system of equations describing the conservation of momentum Eq. (6) [Ha=0], the shear stress Eq. (7), and the energy equation are given as [28, 29]:
$R e\left[\frac{\partial u}{\partial t}+W e^{\alpha} \frac{\partial^{\alpha+1} u}{\partial t^{\alpha+1}}\right]\left[\frac{\partial u}{\partial y} \frac{\partial}{\partial y}\left(\frac{1}{1+\theta_{e} \theta}\right)+\lambda_{r}^{\beta} \frac{\partial^{\beta+1} u}{\partial t^{\beta} \partial y} \frac{\partial}{\partial y}\left(\frac{1}{1+\theta_{e} \theta}\right)+\frac{1}{1+\theta_{e} \theta} \frac{\partial^{2} u}{\partial y^{2}}+\frac{1}{1+\theta_{e} \theta} \frac{\partial^{\beta+2} u}{\partial t^{\beta} \partial y^{2}}\right]=0$ (51)
$R e\left[\frac{\partial \theta}{\partial t}+W e^{\alpha} \frac{\partial^{\alpha+1} \theta}{\partial t^{\alpha+1}}\right]\frac{1}{\operatorname{Pr}}\left[\frac{\partial^{2} \theta}{\partial y^{2}}+W e^{\alpha} \frac{\partial^{\alpha+2} \theta}{\partial t^{\alpha} \partial y^{2}}\right]E c\left[\frac{1}{1+\theta_{e} \theta}\left(\frac{\partial u}{\partial y}\right)^{2}+\frac{\lambda_{r}^{\beta}}{1+\theta_{e} \theta} \frac{\partial^{\beta+1} u}{\partial t^{\beta} \partial y} \frac{\partial u}{\partial y}\right]=0$ (52)
$\left[1+W e^{\alpha} \frac{\partial^{\alpha}}{\partial t^{\alpha}}\right] \tau=\frac{1}{1+\theta_{e} \theta}\left[\frac{\partial u}{\partial y}+\lambda_{r}^{\beta} \frac{\partial^{\beta+1} u}{\partial t^{\beta} \partial y}\right]$ (53)
where, θ_{e}=γΔT, the temperature ratio, $E c=\mathbf{U}_{0}^{2}(c \Delta T)^{1}$, the Eckert number, Pr=cμ_{0} k^{1}, the reference Prandtl number, θ=(TT_{0})ΔT^{1}, the dimensionless temperature are dimensionless parameters. The initial and boundary conditions are given as in Eq. (45) and Eq. (46) and
$\theta(y, 0)=\frac{\partial}{\partial t} \theta(y, 0)=0 \forall y>0$ (54)
$\theta(0, t)=1, \quad \theta(y, t)=0$ as $y \rightarrow \infty, t>0$ (55)
We seek a solution in the form of Eq. (22) for each dependent variable and apply the quasi–linearization method to linearize the system of equations (see Bellman and Kalaba [30], Motsa et al. [31]):
$\begin{aligned}\left[\operatorname{Re}\left[\left(\mathrm{D}_{t}^{1} \otimes I_{y}\right)+W e^{\alpha}\left(D_{t}^{\alpha+1} \otimes I_{y}\right)\right]+a_{0 r}\left(I_{t} \otimes D_{y}^{2}\right)\right.\\+& a_{1 r}\left(I_{t} \otimes \mathrm{D}_{y}\right) \\+a_{2 r}\left(\mathrm{D}_{t}^{\beta} \otimes \mathrm{D}_{y}^{2}\right)+a_{3 r}\left(\mathrm{D}_{t}^{\beta}\right.& \\\left.\left.\otimes \mathrm{D}_{y}^{1}\right)\right] \operatorname{vec}\left(\vec{u}_{r+1}\right) \\+\left[a_{4 r}\left(I_{t} \otimes \mathrm{D}_{y}^{1}\right)\right.& \\\left.+a_{5 r}\left(I_{t} \otimes I_{y}\right)\right] \operatorname{vec}\left(\vec{\theta}_{r+1}\right) &=v e c\left(\vec{R}_{1 r}\right) \end{aligned}$ (56)
$\left[\operatorname{Re}\left(\mathrm{D}_{t}^{1} \otimes I_{y}\right)+\operatorname{Re} W e^{\alpha}\left(\mathrm{D}_{t}^{\alpha+1} \otimes I_{y}\right)\frac{1}{\operatorname{Pr}}\left(I_{t} \otimes \mathrm{D}_{y}^{2}\right)\right.$
$\frac{W e^{\alpha}}{\operatorname{Pr}}\left(\mathrm{D}_{t}^{\alpha} \otimes \mathrm{D}_{y}^{2}\right)+b_{2 r}\left(I_{t}\right.$$\left.\left.\otimes I_{y}\right)\right] \operatorname{vec}\left(\vec{\theta}_{r+1}\right)$
$+\left[b_{0 r}\left(I_{t} \otimes \mathrm{D}_{y}^{1}\right)+b_{1 r}\left(\mathrm{D}_{t}^{\beta}\right.\right.$$\left.\left.\otimes \mathrm{D}_{y}^{1}\right)\right] \operatorname{vec}\left(\vec{u}_{r+1}\right)=\operatorname{vec}\left(\vec{R}_{2 r}\right)$ (57)
Table 2. Residuals for the dependent variables in the problem in Section 4.2
(α,β) 
(N_{t},N_{y}) 
Res_{u} 
Res_{θ} 
Res_{τ} 
(0.25,0.5) 
(7,10) 
2.1518e14 
1.0823e13 
9.8749e16 

(10,15) 
1.7957e13 
6.2928e13 
6.8437e15 

(12,20) 
2.0137e12 
5.9289e12 
3.9240e14 
(0.5,0.75) 
(7,10) 
3.2932e14 
1.7114e13 
1.0634e15 

(10,15) 
4.0798e13 
9.5568e13 
4.3160e15 

(12,20) 
5.1502e12 
8.9247e12 
4.9766e14 
(0.75,1) 
(7,10) 
4.4381e14 
1.3593e13 
2.2664e15 

(10,15) 
1.1795e12 
1.3475e12 
2.0567e14 

(12 ,20) 
1.4751e11 
1.9498e11 
3.6492e14 
(1,1) 
(7,10) 
7.6759e14 
1.3427e13 
3.5037e15 

(10,15) 
2.0559e12 
3.0306e12 
1.7538e14 
(a)
(b)
Figure 4. Velocity and temperature profiles for nonisothermal flow of the special cases of the fractional OldroydB fluid
(a)
(b)
Figure 5. Velocity profiles of the nonisothermal flow for (a) β= 0.5 and (b) α=0.1
Table 2 shows the maximum residual error for 0≤t≤10 for different values of α and β. In Figure 4, the velocity and temperature profiles are illustrated for the special cases of the fractional Oldroyd–B fluid with α=0.25, β=0.5, Re=1.1, Pr=5, Ec=0.2, θe=0.33 and ω=10+π/2. We consider the cases: We=0.2<λr=1.5, Jeffrey fluid (λr=1.2<We=1.5), Maxwell fluid (We=1.5, λr=0), secondgrade fluid (We=0, λr=1.2) and Newtonian fluid. It can be seen that, as in Section 4.1, the secondorder fluid has the smallest amplitude while the Maxwell fluid has the highest oscillation amplitude. However, for the temperature profiles, the Maxwell fluid reaches its thermal boundary layer region quicker than the other fluids, and the secondgrade fluid has the most extensive thermal boundary layer. Figure 4b also show increased temperature close to the wall. The study of Ishfaq et al. [28] attributes this behaviour to the accumulation of energy in the fluid particles in the vicinity of the wall because the viscous effect is profoundly felt in this region. The effects of different orders of derivatives are shown in Figure 5. Unlike the problem in Section 4.1, where velocity profiles are either strictly increasing or decreasing functions of α and β, in this flow, the velocity profiles intersect. We ascribe this behaviour to the coupling effect and non–isothermic nature of the fluid.
In this study, efficient and accurate numerical solutions for unsteady and unidirectional flows of generalized Oldroyd–B fluids have been obtained using the shifted Chebyshev polynomials of the first kind as the basis function. The two problems presented have a varying degree of complexity, which are mainly non–linearity and coupling. The accuracy of the solutions was investigated through the calculation of the residual error norm. The condition number of the discretization matrices was also presented for various values of the orders of derivatives and truncation of the Chebyshev series. The effects of the fractionalorder derivatives in the constitutive equation of the Oldroyd–B fluid on shear stress and velocity are also discussed. We found that small values of the order of the relaxation term can be used to model the flow of viscoelastic fluid with shortterm memory and slow response to shear force. Since many fluid dynamics problems exist on complex geometries, improvement can be made to accommodate problems with complex geometries.
B_{0} the applied magnetic field, N⋅s / C⋅m
c specific heat capacity, J/ K/ kg
Ec Eckert number
e_{ij} stress tensor, N/m^{2}
Ha Hartmann number
k thermal conductivity, W/ m⋅K
Pr reference Prandtl number
Re Reynolds number
U_{0} reference velocity, m/s
We Weissenberg number
Greek symbols
α,β fractional orders
θ_{e} temperature ratio
λ relaxation time, s
λ_{r} retardation time, s/ dimensionless retardation time
μ viscosity, kg/m/s
μ_{0} reference viscosity, kg/m/s
ρ fluid density, kg/m^{3}
σ electrical conductivity, S/m
τ shear stress, N/m^{2}
ω frequency of oscillation, s^{1}
The vectors R_{1r} and R_{2r} in Eq. (56) and Eq. (57) are defined as
$\begin{aligned} R_{1 r}=a_{0 r} \frac{\partial^{2} u_{r}}{\partial y^{2}}+& a_{1 r} \frac{\partial u_{r}}{\partial y}+a_{2 r} \frac{\partial^{\beta+2} u_{r}}{\partial t^{\beta} \partial y^{2}} \\ &+a_{3 r} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y}+a_{4 r} \frac{\partial \theta_{r}}{\partial y}+a_{5 r} \theta_{r} \\ &N u_{r} \end{aligned}$ (58)
where,
$\begin{aligned} N u_{r}=\frac{\partial u_{r}}{\partial y} \frac{\partial}{\partial y}\left(\frac{1}{1+\theta_{e} \theta_{r}}\right)\lambda_{r}^{\beta} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y} \frac{\partial}{\partial y}\left(\frac{1}{1+\theta_{e} \theta_{r}}\right) \\ &\frac{1}{1+\theta_{e} \theta_{r}} \frac{\partial^{2} u_{r}}{\partial y^{2}}\frac{1}{1+\theta_{e} \theta_{r}} \frac{\partial^{\beta+2} u_{r}}{\partial t^{\beta} \partial y^{2}} \end{aligned}$
$a_{0 r}=\left(1+\theta_{e} \theta_{r}\right)^{1}, \quad a_{1 r}=\theta_{e}\left(1+\theta_{e} \theta_{r}\right)^{2} \frac{\partial \theta_{r}}{\partial y}$
$a_{2 r}=\left(1+\theta_{e} \theta_{r}\right)^{1}$$a_{3 r}=\theta_{e} \lambda_{r}^{\beta}\left(1+\theta_{e} \theta_{r}\right)^{2} \frac{\partial \theta_{r}}{\partial y}$
$a_{4 r}=\theta_{e}\left(1+\theta_{e} \theta_{r}\right)^{2}\left(\frac{\partial u_{r}}{\partial y}+\lambda_{r}^{\beta} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y}\right)$
$\begin{aligned} a_{5 r}=2 \theta_{e}^{2}\left(1+\theta_{e} \theta_{r}\right)^{3}\left(\frac{\partial \theta_{r}}{\partial y} \frac{\partial u_{r}}{\partial y}+\lambda_{r}^{\beta} \frac{\partial \theta_{r}}{\partial y} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y}\right) \\ &+\theta_{e}\left(1+\theta_{e} \theta_{r}\right)^{2}\left(\frac{\partial^{2} u_{r}}{\partial y^{2}}+\frac{\partial^{\beta+2} u_{r}}{\partial t^{\beta} \partial y^{2}}\right) \end{aligned}$
and
$R_{2 r}=b_{0 r} \frac{\partial u_{r}}{\partial y}+b_{1 r} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y}+b_{2 r} \theta_{r}N \theta r$ (59)
Here,
$N \theta_{r}=E c\left[\frac{1}{1+\theta_{e} \theta_{r}}\left(\frac{\partial u_{r}}{\partial y}\right)^{2}+\frac{\lambda_{r}^{\beta}}{1+\theta_{e} \theta_{r}} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y} \frac{\partial u_{r}}{\partial y}\right]$
$b_{0 r}=2 E c\left(1+\theta_{e} \theta_{r}\right)^{1} \frac{\partial u_{r}}{\partial y}\lambda_{r}^{\beta}\left(1+\theta_{e} \theta_{r}\right)^{1} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y}$
$b_{1 r}=\lambda_{r}^{\beta}\left(1+\theta_{e} \theta_{r}\right)^{1} \frac{\partial u_{r}}{\partial y}$
$b_{2 r}=\theta_{e}\left(1+\theta_{e} \theta_{r}\right)^{1}\left(E c\left(\frac{\partial u_{r}}{\partial y}\right)^{2}+\lambda_{r}^{\beta} \frac{\partial^{\beta+1} u_{r}}{\partial t^{\beta} \partial y} \frac{\partial u_{r}}{\partial y}\right)$
[1] Oldroyd, J. (1958). NonNewtonian effects in steady motion of some idealized elasticoviscous liquids. Proceedings of the Royal Society London Series A, 245: 278297.
[2] Maxwell, J.C. (2003). On the dynamical theory of gases: The Kinetic Theory of Gases: An anthology of classic papers with historical commentary. World Scientific, 197261.
[3] Tarasov, V.E. (2013). Review of some promising fractional physical models. International Journal of Modern Physics B, 27(9): 1330005. http://dx.doi.org/10.1142/S0217979213300053
[4] Jamil, M., Khan, N., Zafar, A. (2011). Translational flows of an OldroydB fluid with fractional derivatives. Computers & Mathematics with Applications, 62(3): 15401553. http://dx.doi.org/10.1016/j.camwa.2011.03.090
[5] Bazhlekova, E., Bazhlekov, I. (2017). Stokes’ first problem for viscoelastic fluids with a fractional Maxwell model. Fractal and Fractional, 1(1): 7. http://dx.doi.org/10.3390/fractalfract1010007
[6] Le Roux, C. (2014). On flows of viscoelastic fluids of Oldroyd type with wall slip. Journal of Mathematical Fluid Mechanics, 16(2): 335350. http://dx.doi.org/10.1007/s0002101301599
[7] Wang, Q. (2006). Numerical solutions for fractional KdV–Burgers equation by Adomian decomposition method. Applied Mathematics and Computation, 182(2): 10481055. http://dx.doi.org/10.1016/j.amc.2006.05.004
[8] Diethelm, K., Ford, N.J., Freed, A.D. (2002). A predictorcorrector approach for the numerical solution of fractional differential equations. Nonlinear Dynamics, 29: 322.
[9] Erjaee, G., Taghvafard, H., Alnasr M. (2011). Numerical solution of the high thermal loss problem presented by a fractional differential equation. Communications in Nonlinear Science and Numerical Simulation, 169(3): 13561362. http://dx.doi.org/10.1016/j.cnsns.2010.06.031
[10] ElSayed, A., ElKalla, I., Ziada, E. (2010). Analytical and numerical solutions of multiterm nonlinear fractional orders differential equations. Applied Numerical Mathematics, 60(8): 788797. http://dx.doi.org/10.1016/j.apnum.2010.02.007
[11] Momani, S., Odibat, Z. (2007) Comparison between the homotopy perturbation method and the variational iteration method for linear fractional partial differential equations. Computers & Mathematics with Applications, 54(78): 910919. http://dx.doi.org/10.1016/j.camwa.2006.12.037
[12] Momani, S., Odibat, Z. (2007). Homotopy perturbation method for nonlinear partial differential equations of fractional order. Physics Letters A, 365(56): 345350. http://dx.doi.org/10.1016/j.physleta.2007.01.046
[13] Doha, E., A. Bhrawy, A., EzzEldien, S. (2011). Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations. Applied Mathematical Modelling, 35(12): 56625672. http://dx.doi.org/10.1016/j.apm.2011.05.011
[14] Atabakzadeh, M., Akrami, M., Erjaee, G. (2013). Chebyshev operational matrix method for solving multiorder fractional ordinary differential equations. Applied Mathematical Modelling, 37(2021): 89038911. http://dx.doi.org/10.1016/j.apm.2013.04.019
[15] Liu, J., Li, X., Wu, L. (2016). An operational matrix of fractional differentiation of the second kind of Chebyshev polynomial for solving multiterm variable order fractional differential equation. Mathematical Problems in Engineering. http://dx.doi.org/10.1155/2016/7126080
[16] Benattia, M.E., Belghaba, K. (2017). Numerical solution for solving fractional differential equations using shifted Chebyshev wavelet. General Letters in Mathematics, 3: 101110.
[17] Bhrawy, A.H., Alghamdi, M.A. (2013). The operational matrix of Caputo fractional derivatives of modified generalized Laguerre polynomials and its applications. Advances in Difference Equations, 307. http://dx.doi.org/10.1186/168718472013307
[18] Baseri, A., Babolian, E., Abbasbandy, S. (2017). Normalized Bernstein polynomials in solving spacetime fractional diffusion equation. Advances in Difference Equations, 346. http://dx.doi.org/10.1186/s1366201714011
[19] Podlubny, I. (1998). Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Elsevier.
[20] Atangana, A. (2017). Fractional operators with constant and variable order with application to geohydrology. Academic Press.
[21] Abramowitz, M., Stegun, I.A. (1965). Handbook of mathematical functions: with formulas, graphs, and mathematical tables, Courier Corporation.
[22] Ahmadi Darani, M., Saadatmandi, A. (2017). The operational matrix of fractional derivative of the fractionalorder Chebyshev functions and its applications. Computational Methods for Differential Equations, 5: 6787.
[23] Ahmadi Darani, M., Nasiri, M. (2013). A fractional type of the Chebyshev polynomials for approximation of solution of linear fractional differential equations. Computational Methods for Differential Equations, 1: 96107.
[24] Shao, W., Wu, X., Chen, S. (2012). Chebyshev tau meshless method based on the integration–differentiation for biharmonictype equations on irregular domain. Engineering Analysis with Boundary Elements, 36(12): 17871798. http://dx.doi.org/10.1016/j.enganabound.2012.06.005
[25] Bhrawy, A.H. (2016). A highly accurate collocation algorithm for 1+1 and 2+1 fractional percolation equations. Journal of Vibration and Control, 22(9): 22882310. http://dx.doi.org/10.1177/1077546315597815
[26] Jamil, M., Khan, A.N., Shahid, N. (2013). Fractional magnetohydrodynamics OldroydB fluid over an oscillating plate. Thermal Science, 17(4): 9971011. http://dx.doi.org/10.2298/TSCI110731140J
[27] Lai, F., Kulacki, F. (1990). The effect of variable viscosity on convective heat transfer along a vertical surface in a saturated porous medium. International Journal of Heat and Mass Transfer, 33(5): 10281031. http://dx.doi.org/10.1016/00179310(90)900848
[28] Ishfaq, N., Khan, W., Khan, Z. (2017). The Stokes’ second problem for nanofluids. Journal of King Saud UniversityScience. 31(1): 61. http://dx.doi.org/10.1016/j.jksus.2017.05.001
[29] Zhao, J., Zheng, L., Zhang, X., Liu, F. (2016). Unsteady natural convection boundary layer heat transfer of fractional Maxwell viscoelastic fluid over a vertical plate. International Journal of Heat and Mass Transfer, 97: 760766. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.02.059
[30] Bellman, R.E., Kalaba, R.E. (1965). Quasilinearization and nonlinear boundaryvalue problems.
[31] Motsa, S.S., Magagula, V., Sibanda, P. (2014). A bivariate Chebyshev spectral collocation quasilinearization method for nonlinear evolution parabolic equations. The Scientific World Journal. http://dx.doi.org/10.1155/2014/581987