MHD Laminar Boundary Layer Flow of a Jeffrey Fluid Past a Vertical Plate Influenced by Viscous Dissipation and a Heat Source/Sink

Abstract

This study investigates the effects of viscous dissipation and a heat source or sink on the magneto-hydrodynamic laminar boundary layer flow of a Jeffrey fluid past a vertical plate. The governing boundary layer non-linear partial differential equations are reduced to non-linear ordinary differential equations using suitable similarity transformations. The resulting system of dimensionless differential equations is then solved numerically using the bivariate spectral quasi-linearisation method. The effects of some physical parameters that include the Schmidt number, Eckert number, radiation parameter, magnetic field parameter, heat generation parameter, and the ratio of relaxation to retardation times on the velocity, temperature, and concentration profiles are presented graphically. Additionally, the influence of some physical parameters on the skin friction coefficient, local Nusselt number, and the local Sherwood number are displayed in tabular form.

1. Introduction

The study of fluids, either stationary or in motion has vast applications in science, biology, physiology, medicine, and engineering. Basically, fluids can be grouped into two classes, Newtonian and non-Newtonian fluids. Newtonian fluids are those fluids that obey Newton’s law of viscosity. Those fluids that deviate from the law are said to be non-Newtonian. Typical non-Newtonian fluids include blood, saliva, starch solutions, toothpaste, and sauces. Due to the fact that non-Newtonian fluids cannot be described by a single constitutive equation like the Navier–Stokes equations, they tend to be more complicated than Newtonian fluids. Although complicated, non-Newtonian fluids have found more industrial applications than Newtonian fluids. Recently, there is much attention which has been paid on the study of non-Newtonian fluid flow. The non-Newtonian fluids can be characterised as being of the differential, rate or integral type, Cioranescu et al. [1]. Some of the non-Newtonian models which have been studied include the Jeffrey model [2,3], the Maxwell model [4,5], the Oldroyd-B model [6,7], and the Herschel–Bulkley model [8,9].

One non-Newtonian fluid model of interest in this study is the Jeffrey model which is of the rate type. The Jeffrey fluid expresses the influence of the ratio of relaxation to retardation times and its constitutive equation can be reduced to that of a Newtonian fluid as a special case. Unlike viscous fluid models, the Jeffrey model can be used in non-Newtonian fluids to describe the stress relaxation property, Kahshan et al. [10]. Due to its viscoelastic properties, the Jeffrey fluid has found some industrial applications in polymer productions. The model has recently attracted the attention of many researchers since it gives better approximations to most physiological fluids. The Jeffrey fluid model was successfully used to model the peristaltic flow of chyme in the small intestines by Akbar et al. [11]. Additionally, Sharma et al. [12] used the Jeffrey fluid to model the flow of blood in narrow arteries. Nallapu and Radhakrishnamacharya [13] did an investigation of the effects of a magnetic field in the flow of the Jeffrey fluid in narrow tubes in a porous medium. Some other researches which have been completed on the Jeffrey model include the works by: Ellahi et al. [14], Khan et al. [15], Vaidya et al. [16].

There is quite a substantial amount of research which has been performed on the MHD boundary layer flow of a Jeffrey fluid. A study on magneto-hydrodynamic boundary layer flow of a Jeffrey fluid over a sheet in a porous medium was conducted by Ahmad and Ishak [17]. Nadeem et al. [18] did a numerical study of the laminar boundary layer flow of a Jeffrey fluid in which they investigated the effects of thermal radiation when the fluid is flowing over a surface that stretches exponentially. The homotopy analysis approach was employed by Hayat et al. [19] to study the influence of some embedded parameters on the laminar boundary layer flow of a Jeffrey fluid. Tlili [20] investigated the heat generation and MHD effects on an incompressible two-dimensional laminar boundary layer Jeffrey fluid flow.

Some significant work has been completed on the non-Newtonian fluid flow past a vertical plate. This is attributed to its wide applications in engineering, technology, as well as in the manufacturing industry. Babu et al. [21] considered the influence of the Hall currents in a chemically reactive multivariate Jeffrey fluid which is flowing through a vertical plate. A study on the effects of viscous dissipation on an MHD Jeffrey fluid flowing through a vertical channel was completed by Selvi and Muthuraj [22]. Amanulla et al. [23] discussed the effects of heat transfer in a Jeffrey fluid from an inclined vertical plate. An investigation on the flow of an incompressible non-Newtonian Jeffrey fluid past a vertical porous plate was completed by Prasad et al. [24] using an implicit finite difference method.

The study of flow of non-Newtonian fluids influenced by a heat source or sink has found applications in the disposal of radioactive waste material, food stuffs storage, removal of heat from nuclear fuel debris, and many others, Hayat et al. [25]. Nisar et al. [26] did an analysis of the flow of an electrically conducting fluid flowing over a stretching surface with a uniform heat source. A Cattaneo–Christov model with double stratification, thermal relaxation, and a heat source was used to investigate heat and mass transfer in a Jeffrey fluid by Hussain et al. [27]. Qasim [28] investigated both heat and mass transfer in a Jeffrey fluid flowing over a stretching sheet with a heat source.

In this current work, a numerical technique called the bivariate spectral quasi-linearisation method (BSQLM), Motsa et al. [29], was used to perform a mathematical investigation of the laminar boundary layer flow of a Jeffrey fluid. The non-Newtonian fluid flows past a vertical plate with the flow being influenced by viscous dissipation and a heat source or sink. We intend to implement the BSQLM in this current study as it is reported by Magagula et al. [30] that the method is computationally efficient and gives solutions that are more uniformly accurate than traditional methods like the finite difference methods. Although the BSQLM does not handle periodic boundary conditions properly, the method is considered because it is easy to implement and it takes little time to converge to the exact solution with few grid points, Rai and Mondal [31]. There is quite a number of studies that have been completed in analysing fluid flow problems using the BSQLM. Goqo et al. [32] used the BSQLM to obtain the convergent solutions of the model equations that describe the viscous nanofluid flowing over a porous wedge. Additionally, Motsa and Ansari et al. [33] used the BSQLM to investigate unsteady boundary layer flow of an Oldroyd-B nanofluid towards a stretching sheet with variable thermal conductivity.

2. Problem Statement and Mathematical Formulation

In this study, we consider the magneto-hydrodynamic laminar boundary layer flow of a Jeffrey fluid past a vertical porous plate, with the flow being influenced by viscous dissipation and a heat source or sink. The x-axis is measured along the plate in the upward direction and the y-axis is measured normal to the plate in the outward direction. Initially, both the surface of the plate and the Jeffrey fluid are at rest at a uniform temperature $T_w$ and concentration $C_w$, respectively. The respective ambient temperature and concentration of the fluid far way from the plate are given by $T_{\infty }$ and $C_{\infty }$. The acceleration due to gravity g acts vertically downwards. The geometry of the problem is shown in Figure 1.

The Jeffrey fluid constitutive equation is given by [34,35] as,

$$\begin{array}{ccc}& & \mathbf{T}=-p\mathbf{I}+\mathbf{S},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & \mathbf{S}=\frac{\mu }{1+\lambda }\left[{\mathbf{A}}_1+{\lambda }_1\frac{d}{dt}{\mathbf{A}}_1\right],\hfill \end{array}$$

(2)

where $\mathbf{S}$ is the extra-stress tensor, ${\mathbf{A}}_1$ is the first Rivlin–Erickson tensor, $\lambda$ is the ratio of relaxation to retardation times, ${\lambda }_1$ is the retardation time and $\mu$ is the viscosity. After incorporating the Jeffrey fluid stress tensor in the momentum equation, the model is given by, Abdul Gaffar et al. [36]

$$\begin{array}{ccc}& & \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{\nu }{1+\lambda }\left(\frac{{\partial }^{2}u}{\partial y^2}\right)+g{\beta }_T(T-T_{\infty })+g{\beta }_C(C-C_{\infty })-\frac{\sigma B_0^2}{\rho }u\hfill \\ & & +\frac{\nu {\lambda }_1}{1+\lambda }\left(u\frac{{\partial }^{3}u}{\partial x\partial y^2}+\frac{\partial u}{\partial y}\frac{{\partial }^{2}u}{\partial x\partial y}-\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial y^2}+v\frac{{\partial }^{3}u}{\partial y^3}\right),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2}\pm \frac{Q}{\rho c_p}(T-T_{\infty })+\frac{\mu }{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^2+\frac{\sigma B_0^2}{\rho c_p}{u}^2-\frac{1}{\rho c_p}\frac{\partial q_r}{\partial y},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_m\frac{{\partial }^{2}C}{\partial y^2}-k_r(C-C_{\infty }),\hfill \end{array}$$

(6)

where the velocities in the x and y directions are u and v, respectively. T is the fluid temperature, C is the fluid concentration, $D_m$ is mass diffusivity, $k_r$ is the chemical reaction rate constant, $\nu$ is kinematic viscosity, $q_r$ is the radiative heat flux, $\rho$ is fluid density. Additionally, g, ${\beta }_T$, ${\beta }_C$ and $\alpha$ are the acceleration due to gravity, thermal expansion coefficient, solutal expansion coefficient, and thermal diffusivity, respectively. The boundary conditions are chosen in such a way that there is no flow across the plate and due to viscosity, there is no slip condition at the plate. These boundary conditions for the temperature, concentration, and velocity are given by:

$$\begin{array}{ccc}& & u=0,v=0,T=T_w,C=C_{w}aty=0,\hfill \\ \hfill & & u\to 0,\frac{\partial u}{\partial y}\to 0,C\to C_{\infty },T\to T_{\infty }asy\to \infty .\hfill \end{array}$$

(7)

Under the Rosseland approximation, the radiative heat flux $q_r$ has the form, El-Aziz and Yahya [37]

$$\begin{array}{c}\hfill q_r=-\frac{4{\sigma }^{*}}{3k^{*}}\frac{\partial T^4}{\partial y},\end{array}$$

(8)

where ${\sigma }^{*}$ is the Stefan–Boltzmann constant and $k^{*}$ is the mean absorption coefficient. Assuming that the temperature differences within the Jeffrey fluid flow are sufficiently small, Taylor’s series expansion of $T^4$ about $T_{\infty }$ gives

$$\begin{array}{c}\hfill T^4\simeq 4T_{\infty }^{3}T-3T_{\infty }^4.\end{array}$$

(9)

From Equation (8), the derivative of $q_r$ with respect to y is

$$\begin{array}{c}\hfill \frac{\partial q_r}{\partial y}=-\frac{16{\sigma }^{*}{T}_{\infty }^3}{3k^{*}}\frac{{\partial }^{2}T}{\partial y^2},\end{array}$$

(10)

which reduces Equation (5) to the form

$$\begin{array}{c}\hfill u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2}\pm \frac{Q}{\rho c_p}(T-T_{\infty })+\frac{\mu }{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^2+\frac{\sigma B_0^2}{\rho c_p}{u}^2+\frac{16{\sigma }^{*}{T}_{\infty }^3}{3\rho c_p{k}^{*}}\frac{{\partial }^{2}T}{\partial y^2}.\end{array}$$

(11)

Similarity Transformations

The transformation of the system of Equations (4)–(6) into a dimensionless system is achieved by the application of the similarity transformations, Abdul Gaffar et al. [36]:

$$\begin{array}{ccc}& & \xi ={\left(\frac{x}{b}\right)}^{\frac{1}{2}},\eta =ay{x}^{\frac{1}{4}},a={\left(\frac{g{\beta }_t(T_w-T_{\infty })}{4{\nu }^2}\right)}^{\frac{1}{4}},\psi =4\nu a{x}^{\frac{3}{4}}f(\eta ,\xi ),\hfill \\ & & \theta (\eta ,\xi )=\frac{T-T_{\infty }}{T_w-T_{\infty }},\varphi (\eta ,\xi )=\frac{C-C_{\infty }}{C_w-C_{\infty }}.\hfill \end{array}$$

(12)

The similarity transformations defined in Equation (12) when applied to the momentum, energy, and concentration equations yield

$$\begin{array}{ccc}& & {\mathsf{\Lambda }}_f[f,\theta ,\varphi ]:=\frac{f^{‴}}{1+\lambda }+3f{f}^{″}-2{(f^{\prime })}^2+2\xi (f^{″}{f}_{\xi }-f^{\prime }{f}_{\xi }^{\prime })+\theta +K\varphi -M{f}^{\prime }\hfill \\ & & +2\xi \frac{De}{1+\lambda }(f^{\prime }{f}_{\xi }^{‴}+f^{″}{f}_{\xi }^{″}-f^{‴}{f}_{\xi }^{\prime }-f^{iv}{f}_{\xi })+\frac{De}{1+\lambda }(f^{″}{}^2-2f^{\prime }{f}^{‴}-3f{f}^{iv})=0,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & {\mathsf{\Lambda }}_{\theta }[f,\theta ,\varphi ]:=(1+\frac{4}{3}{R}_d){\theta }^{″}+3Prf{\theta }^{\prime }+Pr\Delta \theta +PrEc({(f^{″})}^2+M{(f^{\prime })}^2)+2\xi Pr({\theta }^{\prime }{f}_{\xi }-f^{\prime }{\theta }_{\xi })=0,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\mathsf{\Lambda }}_{\varphi }[f,\theta ,\varphi ]:={\varphi }^{″}+3Scf{\varphi }^{\prime }-\gamma Sc\varphi -2\xi Sc(f^{\prime }{\varphi }_{\xi }-{\varphi }^{\prime }{f}_{\xi })=0,\hfill \end{array}$$

(15)

where $\gamma =\frac{k_r{x}^{\frac{1}{2}}}{\nu B^2}$ is the chemical reaction parameter, $R_d=\frac{\nu k^{*}}{4T_{\infty }^3{\sigma }^{*}}$ is the radiation parameter, $K=\frac{{\beta }_C(C_w-C_{\infty })}{{\beta }_T(T_w-T_{\infty })}$ is the concentration to thermal buoyancy ratio, $Pr=\frac{\nu }{\alpha }$ is the Prandtl number, $Sc=\frac{\nu }{D_m}$ is the Schmidt number, $M=\frac{\sigma B_0^2{x}^{\frac{1}{2}}}{\nu C_1^2\rho }$ is the magnetic field parameter, $Ec=\frac{16a^4{\nu }^{2}x}{c_p(T_w-T_{\infty })}$ is the Eckert number, $\Delta =\frac{Q{x}^{\frac{1}{2}}}{\rho \nu a^2{c}_p}$ is the heat generation or absorption parameter and $De=\frac{{\lambda }_1\nu a^2}{x^{1/2}}$ is the Deborah number. The dimensionless boundary conditions take the form

$$\begin{array}{ccc}& & f(0,\xi )=0,f^{\prime }(0,\xi )=0,\theta (0,\xi )=1,\varphi (0,\xi )=1,\hfill \\ & & f^{\prime }(\infty ,\xi )\to 0,f^{″}(\infty ,\xi )\to 0,\theta (\infty ,\xi )\to 0,\varphi (\infty ,\xi )\to 0.\hfill \end{array}$$

(16)

3. Method of Solution

The numerical solutions of the dimensionless system of Equations (13)–(15) subject to boundary conditions (16) are obtained using the BSQLM. The method basically combines the quasi-linearisation method (QLM) and the Chebyshev spectral collocation method (CSCM). The application of a Newton–Raphson based QLM, Alharbey et al. [38], transforms the coupled non-linear differential Equations (13)–(15) into a system of linear iterative differential equations:

$$\begin{array}{ccc}& & a_{0,m}{f}_{m+1}^{iv}+a_{1,m}{f}_{m+1}^{‴}+a_{2,m}{f}_{m+1}^{″}+a_{3,m}{f}_{m+1}^{\prime }+a_{4,m}{f}_{m+1}+a_{5,m}{\theta }_{m+1}+a_{6,m}{\varphi }_{m+1}+a_{7,m}\frac{{\partial }^4{f}_{m+1}}{\partial {\eta }^3\partial \xi }\hfill \\ & & +a_{8,m}\frac{{\partial }^3{f}_{m+1}}{\partial {\eta }^2\partial \xi }+a_{9,m}\frac{{\partial }^2{f}_{m+1}}{\partial \eta \partial \xi }+a_{10,m}\frac{\partial f_{m+1}}{\partial \xi }={\mathcal{R}}_{1,m},\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & b_{0,m}{\theta }_{m+1}^{″}+b_{1,m}{\theta }_{m+1}^{\prime }+b_{2,m}{\theta }_{m+1}+b_{3,m}{f}_{m+1}^{″}+b_{4,m}{f}_{m+1}^{\prime }+b_{5,m}{f}_{m+1}+b_{6,m}\frac{\partial f_{m+1}}{\partial \xi }+b_{7,m}\frac{\partial {\theta }_{m+1}}{\partial \xi }\hfill \\ & & ={\mathcal{R}}_{2,m},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}& & c_{0,m}{\varphi }_{m+1}^{″}+c_{1,m}{\varphi }_{m+1}^{\prime }+c_{2,m}{\varphi }_{m+1}+c_{3,m}{f}_{m+1}^{\prime }+c_{4,m}{f}_{m+1}+c_{5,m}\frac{\partial f_{m+1}}{\partial \xi }+c_{6,m}\frac{\partial {\varphi }_{m+1}}{\partial \xi }={\mathcal{R}}_{3,m}.\hfill \end{array}$$

(19)

The subscripts m and $m+1$ denote the previous and current iteration levels, respectively. The variable coefficients evaluated at the previous iteration level m are expressed as follows:

$$\begin{array}{ccc}& & a_{0,m}=-2\xi \frac{De}{1+\lambda }\frac{\partial f_m}{\partial \xi }-3\frac{De}{1+\lambda }{f}_m,a_{1,m}=\frac{1}{1+\lambda }-2\xi \frac{De}{1+\lambda }\frac{\partial f_m^{\prime }}{\partial \xi }-2\frac{De}{1+\lambda }{f}_m^{\prime },\hfill \\ & & a_{2,m}=3f_m+2\xi \frac{\partial f_m}{\partial \xi }+2\xi \frac{De}{1+\lambda }\frac{\partial f_m^{″}}{\partial \xi }+2\frac{De}{1+\lambda }{f}_s^{″},\hfill \\ & & a_{3,m}=-M-4f_m^{\prime }-2\frac{\partial f_m^{\prime }}{\partial \xi }+2\xi \frac{De}{1+\lambda }\frac{\partial f_m^{‴}}{\partial \xi }-2\frac{De}{1+\lambda }{f}_s^{‴},a_{4,m}=3f_s^{″}-3\frac{De}{1+\lambda }{f}_m^{iv},\hfill \\ & & a_{5,m}=1,a_{6,m}=K,a_{7,m}=2\xi \frac{De}{1+\lambda }{f}_m^{\prime },a_{8,m}=2\xi \frac{De}{1+\lambda }{f}_m^{″},a_{9,m}=-2\xi \frac{De}{1+\lambda }{f}_s^{‴}-2\xi f_m^{\prime },\hfill \\ & & a_{10,m}=2\xi f_m^{″}-2\xi \frac{De}{1+\lambda }{f}_m^{iv},b_{0,m}=1+\frac{4}{3}{R}_d,b_{1,m}=Pr\left(3f_m+2\xi \frac{\partial f}{\partial \xi }\right),b_{2,m}=\Delta Pr,\hfill \\ & & b_{3,m}=2PrEc{f}_m^{″},b_{4,m}=2Pr\left(EcM{f}_m^{\prime }-\xi \frac{\partial {\theta }_m}{\partial \xi }\right),b_{5,m}=3Pr{\theta }_m^{\prime },b_{6,m}=2\xi Pr{\theta }_m^{\prime },b_{7,m}=-2\xi Pr{f}_m^{\prime },\hfill \\ & & c_{0,m}=1,c_{1,m}=Sc\left(3f_m+2\xi \frac{\partial f_m}{\partial \xi }\right),c_{2,m}=-\gamma Sc,c_{3,m}=-2\xi Sc\frac{\partial {\varphi }_m}{\partial \xi },c_{4,m}=3Sc{\varphi }_m^{\prime },\hfill \\ & & c_{5,m}=2\xi Sc{\varphi }_m^{\prime },c_{6,m}=-2\xi Sc{f}_m^{\prime }.\hfill \end{array}$$

Additionally, the right hand side terms, evaluated at the current iteration level, are given by the following definitions

$$\begin{array}{ccc}& & {\mathcal{R}}_{1,m}=3f_m{f}_m^{″}-2{(f_m^{\prime })}^2+2\xi f_m^{″}\frac{\partial f_m}{\partial \xi }-f_m^{\prime }\frac{\partial f_m^{\prime }}{\partial \xi }\hfill \\ & & +\frac{De}{1+\lambda }\left({(f_m^{″})}^2-2f_m^{\prime }{f}_m^{‴}-3f_m{f}_m^{iv}+2\xi \left(f_m^{\prime }\frac{\partial f_m^{‴}}{\partial \xi }+f_m^{″}\frac{\partial f_m^{″}}{\partial \xi }-f_m^{‴}\frac{\partial f_m^{\prime }}{\partial \xi }-f_m^{iv}\frac{\partial f_m}{\partial \xi }\right)\right),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & {\mathcal{R}}_{2,m}=3Pr{f}_m{\theta }_m^{\prime }+PrEcM{(f_m^{\prime })}^2+2Pr\xi {\theta }_m^{\prime }\frac{\partial f_m}{\partial \xi }-2\xi Pr{f}_m^{\prime }\frac{\partial {\theta }_m}{\partial \xi }+PrEc{(f_m^{″})}^2,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & {\mathcal{R}}_{3,m}=3Sc{f}_m{\varphi }_m^{\prime }-2\xi Sc{f}_m^{\prime }\frac{\partial {\varphi }_m}{\partial \xi }+2\xi Sc{\varphi }_m^{\prime }\frac{\partial f_m}{\partial \xi }.\hfill \end{array}$$

(22)

The boundary conditions evaluated at the current iteration level $m+1$ are:

$$\begin{array}{ccc}& & f_{m+1}(0,\xi )=0,f_{m+1}^{\prime }(0,\xi )=0,{\theta }_{m+1}(0,\xi )=1,{\varphi }_{m+1}(0,\xi )=1,\hfill \\ & & f_{m+1}^{\prime }(\infty ,\xi )\to 0,f_{m+1}^{″}(\infty ,\xi )\to 0,{\theta }_{m+1}(\infty ,\xi )\to 0,{\varphi }_{m+1}(\infty ,\xi )\to 0.\hfill \end{array}$$

In order for us to implement the CSCM on the system of Equations (17)–(19), we move the physical domain $[0,L_{\eta }]\times [0,L_{\xi }]$ on the $\eta -\xi$ axis on which Equations (13)–(15) are defined to the computational domain $[-1,1]\times [-1,1]$ on the $x-y$ axis using linear transformations

$$\begin{array}{c}\hfill \left\{\eta =\frac{L_{\eta }}{2}(1+x);\xi =\frac{L_{\xi }}{2}(1+y)\right\}.\end{array}$$

$L_{\eta },L_{\xi }\in {\mathbb{Z}}^{+}$ truncate the semi-infinite physical domains in the $\eta$ and $\xi$ directions, respectively, such that the infinity boundary conditions hold. Approximations of $f_{m+1}$, ${\theta }_{m+1}$ and ${\varphi }_{m+1}$ have the following respective forms

$$\begin{array}{c}\hfill f_{m+1}(\eta ,\xi )=\sum _{j=0}^N\sum _{i=0}^{N^{*}}{f}_{m+1}({\eta }_j,{\xi }_i)L_j(\eta )L_i(\xi ),\\ \hfill {\theta }_{m+1}(\eta ,\xi )=\sum _{j=0}^N\sum _{i=0}^{N^{*}}{\theta }_{m+1}({\eta }_j,{\xi }_i)L_j(\eta )L_i(\xi ),\\ \hfill {\varphi }_{m+1}(\eta ,\xi )=\sum _{j=0}^N\sum _{i=0}^{N^{*}}{\varphi }_{m+1}({\eta }_j,{\xi }_i)L_j(\eta )L_i(\xi ),\end{array}$$

(23)

where the Lagrange interpolating polynomials are defined as:

$$L_j(\eta )=\prod _{\begin{array}{c}k=0\\ k\ne j\end{array}}^N\frac{\eta -{\eta }_k}{{\eta }_j-{\eta }_k},L_i(\xi )=\prod _{\begin{array}{c}k=0\\ k\ne i\end{array}}^{N^{*}}\frac{\eta -{\eta }_k}{{\eta }_i-{\eta }_k},\mathrm{such}\mathrm{that}{L}_j(m_k)=\left\{\begin{array}{c}1,ifj=k,\hfill \\ 0,otherwise.\hfill \end{array}\right.$$

Evaluating the differential Equations (17)–(19) at the Chebyshev–Gauss–Lobatto collocation points ${\displaystyle x_j=cos\left(\frac{\pi j}{N}\right){|}_{j=0}^N}$ and ${\displaystyle y_i=cos\left(\frac{\pi i}{N^{*}}\right){|}_{i=0}^{N^{*}}}$ and Chebyshev differentiation gives

$$\begin{array}{c}({\mathbf{a}}_{0,m,i}{\mathcal{D}}^4+{\mathbf{a}}_{1,m,i}{\mathcal{D}}^3+{\mathbf{a}}_{2,m,i}{\mathcal{D}}^2+{\mathbf{a}}_{3,m,i}\mathcal{D}+{\mathbf{a}}_{4,m,i}){\mathbf{F}}_{m+1,i}+{\mathbf{a}}_{5,m,i}{\theta }_{m+1,i}+{\mathbf{a}}_{6,m,i}{\varphi }_{m+1,i}\hfill \\ +\sum _{k=0}^{N^{*}}({\mathbf{a}}_{7,m,i}{\mathcal{D}}^3+{\mathbf{a}}_{8,m,i}{\mathcal{D}}^2+{\mathbf{a}}_{9,m,i}\mathcal{D}+{\mathbf{a}}_{10,m,i})d_{ik}{\mathbf{F}}_{m+1,k}={\mathbf{R}}_{1,m,i},\hfill \end{array}$$

(24)

$$\begin{array}{c}({\mathbf{b}}_{0,m,i}{\mathcal{D}}^2+{\mathbf{b}}_{1,m,i}\mathcal{D}+{\mathbf{b}}_{2,m,i}){\mathsf{\Theta }}_{m+1,i}+({\mathbf{b}}_{3,m,i}{\mathcal{D}}^2+{\mathbf{b}}_{4,m,i}\mathcal{D}+{\mathbf{b}}_{5,m,i}){\mathbf{F}}_{m+1,i}\hfill \\ +{\mathbf{b}}_{6,m,i}\sum _{k=0}^{N^{*}}{d}_{ik}{\mathbf{F}}_{m+1,k}+{\mathbf{b}}_{7,m,i}\sum _{k=0}^{N^{*}}{d}_{ik}{\theta }_{m+1,k}={\mathbf{R}}_{2,m,i},\hfill \end{array}$$

(25)

$$\begin{array}{c}({\mathbf{c}}_{0,m,i}{\mathcal{D}}^2+{\mathbf{c}}_{1,m,i}\mathcal{D}+{\mathbf{c}}_{2,m,i}){\mathsf{\Phi }}_{m+1,i}+({\mathbf{c}}_{3,m,i}\mathcal{D}+{\mathbf{c}}_{4,m,i}){\mathbf{F}}_{m+1,i}+{\mathbf{c}}_{5,m,i}\sum _{k=0}^{N^{*}}{d}_{ik}{\mathbf{F}}_{m+1,k}\hfill \\ +{\mathbf{c}}_{6,m,i}\sum _{k=0}^{N^{*}}{d}_{iq}{\varphi }_{m+1,k}={\mathbf{R}}_{3,m,i},\hfill \end{array}$$

(26)

where ${\displaystyle \mathcal{D}=\left(\frac{2}{L_{\eta }}\right)\mathbf{D}}$ and ${\displaystyle d=\left(\frac{2}{L_{\xi }}\right)\mathbf{d}}$. $\mathbf{D}$ and $\mathbf{d}$ are the Chebyshev differentiation matrices, Trefethen [39]. The diagonal matrices

$$\begin{array}{c}\hfill {\mathbf{a}}_{k,m,i}=\mathrm{diag}[a_{k,m}({\eta }_0,{\xi }_i),a_{k,m}({\eta }_1,{\xi }_i),\cdots ,a_{k,m}({\eta }_{N-1},{\xi }_i),a_{k,m}({\eta }_N,{\xi }_i)],k=0,1,\cdots ,10,\\ \hfill {\mathbf{b}}_{k,m,i}=\mathrm{diag}[b_{k,m}({\eta }_0,{\xi }_i),b_{k,m}({\eta }_1,{\xi }_i),\cdots ,b_{k,m}({\eta }_{N-1},{\xi }_i),b_{k,m}({\eta }_N,{\xi }_i)],k=0,1,\cdots ,7,\\ \hfill {\mathbf{c}}_{k,m,i}=\mathrm{diag}[c_{k,m}({\eta }_0,{\xi }_i),c_{k,m}({\eta }_1,{\xi }_i),\cdots ,c_{k,m}({\eta }_{N-1},{\xi }_i),c_{k,m}({\eta }_N,{\xi }_i)],k=0,1,\cdots ,6.\end{array}$$

In matrix-vector form, the system of Equations (24)–(26) can be written as

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}M{11}_{i,j}& M{12}_{i,j}& M{13}_{i,j}\\ M{21}_{i,j}& M{22}_{i,j}& M{23}_{i,j}\\ M{31}_{i,j}& M{32}_{i,j}& M{33}_{i,j}\end{array}\right]\left[\begin{array}{c}{\mathbf{F}}_{m+1,i}\\ {\mathsf{\Theta }}_{m+1,i}\\ {\mathsf{\Phi }}_{m+1,i}\end{array}\right]=\left[\begin{array}{c}{\mathbf{R}}_{1,m,i}\\ {\mathbf{R}}_{2,m,i}\\ {\mathbf{R}}_{3,m,i}\end{array}\right],\end{array}$$

(27)

where

$$\begin{array}{c}\hfill {\mathbf{F}}_{m+1,i}={\left[\begin{array}{c}{f}_{m+1}({\eta }_0,{\xi }_i),f_{m+1}({\eta }_1,{\xi }_i),\cdots ,f_{m+1}({\eta }_{N-1},{\xi }_i),f_{m+1}({\eta }_N,{\xi }_i)\end{array}\right]}^T,\\ \hfill {\mathbf{R}}_{1,m,i}={\left[\begin{array}{c}{\mathcal{R}}_{1,m}({\eta }_0,{\xi }_i),{\mathcal{R}}_{1,m}({\eta }_1,{\xi }_i),\cdots ,{\mathcal{R}}_{1,m}({\eta }_{N-1},{\xi }_i),{\mathcal{R}}_{1,m}({\eta }_N,{\xi }_i)\end{array}\right]}^T.\end{array}$$

${\mathsf{\Theta }}_{m+1,i}$ and ${\mathsf{\Phi }}_{m+1,i}$ have a similar form as ${\mathbf{F}}_{m+1,i}$, whilst ${\mathbf{R}}_{2,m,i}$ and ${\mathbf{R}}_{3,m,i}$ take the same form as ${\mathbf{R}}_{1,m,i}$. When $i=j$,

$$\begin{array}{ccc}& & M{11}_{i,j}={\mathbf{a}}_{0,m,i}{\mathbf{D}}^{\mathbf{4}}+{\mathbf{a}}_{1,m,i}{\mathbf{D}}^{\mathbf{3}}+{\mathbf{a}}_{2,m,i}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{a}}_{3,m,i}\mathbf{D}+{\mathbf{a}}_{4,m,i}+d_{ii}({\mathbf{a}}_{7,m,i}{\mathbf{D}}^{\mathbf{3}}+{\mathbf{a}}_{8,m,i}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{a}}_{9,m,i}\mathbf{D}+{\mathbf{a}}_{10,m,i})\hfill \\ & & M{12}_{i,j}=a_{5,m,i}\mathbf{I},M{13}_{i,j}=a_{6,m,i}\mathbf{I},M{21}_{i,j}={\mathbf{b}}_{3,m,i}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{b}}_{4,m,i}\mathbf{D}+{\mathbf{b}}_{5,m,i}+d_{ii}{\mathbf{b}}_{6,m,i},\hfill \\ & & M{22}_{i,j}={\mathbf{b}}_{0,m,i}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{b}}_{1,m,i}\mathbf{D}+b_{2,m,i}\mathbf{I}+d_{ii}{\mathbf{b}}_{7,m,i},M_{23}=\mathbf{0},M{31}_{i,j}={\mathbf{c}}_{3,m,i}\mathbf{D}+{\mathbf{c}}_{4,m,i}+d_{ii}{\mathbf{c}}_{5,m,i},\hfill \\ & & M{32}_{i,j}=\mathbf{0},M{33}_{i,j}={\mathbf{c}}_{0,m,}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{c}}_{1,m,}\mathbf{D}+c_{2,m,i}\mathbf{I}+d_{ii}{\mathbf{c}}_{6,m,i},\hfill \end{array}$$

and when $i\ne j$

$$\begin{array}{ccc}& & M{11}_{i,j}=d_{ij}({\mathbf{a}}_{7,m,i}{\mathbf{D}}^{\mathbf{3}}+{\mathbf{a}}_{8,m,i}{\mathbf{D}}^{\mathbf{2}}+{\mathbf{a}}_{9,m,i}\mathbf{D}+{\mathbf{a}}_{10,m,i}),M{12}_{i,j}=M{13}_{i,j}=\mathbf{0},M{21}_{i,j}=d_{ij}{\mathbf{b}}_{6,m,i},\hfill \\ & & M{22}_{i,j}=d_{ij}{\mathbf{b}}_{7,m,i},M{23}_{i,j}=\mathbf{0},M{31}_{i,j}=d_{ij}{\mathbf{c}}_{5,m,i},M{32}_{i,j}=\mathbf{0},M{33}_{i,j}=d_{ij}{\mathbf{c}}_{6,m,i}.\hfill \end{array}$$

The engineering design quantities of physical interest to be discussed in this work are the skin-friction coefficient $(C_f)$ which is a measure of shear stress at the plate, the Nusselt number $(Nu)$ which is a measure of the rate of heat transfer and the Sherwood number $(Sh)$ which measures the rate of mass transfer at the plate.

4. Results and Discussions

Presented in this section are the numerical results of the BSQLM algorithm for solving MHD laminar boundary layer flow of a Jeffrey fluid past a vertical plate. All the numerical simulations were completed using MATLAB 2020. Unless otherwise stated, the default parameters considered in this work are: $Sc=0.6$, $Pr=0.71$, $\lambda =0.2$, $\Delta =0.1$, $\gamma =0.1$, $K=0.5$, $R_d=0.5$, $E_c=0.1$, and $M=0.1$. The method was tested for both convergence and accuracy using the infinity norms of both the absolute and residual errors, respectively. The formulae for the convergence error infinity norms are defined as:

$$\begin{array}{c}\hfill \mathrm{E}\left[\mathrm{f}\right]=\left|\right|{\mathsf{\Lambda }}_f[f_{m+1},{\theta }_{m+1},{\varphi }_{m+1}]-{\mathsf{\Lambda }}_f[f_m,{\theta }_m,{\varphi }_m]{\left|\right|}_{\infty },\\ \hfill \mathrm{E}\left[\theta \right]=\left|\right|{\mathsf{\Lambda }}_{\theta }[f_{m+1},{\theta }_{m+1},{\varphi }_{m+1}]-{\mathsf{\Lambda }}_{\theta }[f_m,{\theta }_m,{\varphi }_m]{\left|\right|}_{\infty },\\ \hfill \mathrm{E}\left[\varphi \right]=\left|\right|{\mathsf{\Lambda }}_{\varphi }[f_{m+1},{\theta }_{m+1},{\varphi }_{m+1}]-{\mathsf{\Lambda }}_{\varphi }[f_m,{\theta }_m,{\varphi }_m]{\left|\right|}_{\infty }.\end{array}$$

(28)

Figure 2 displays the convergence error norms plotted against the number of iterations. We can infer that the error norms decrease with increasing number of iterations. This confirms convergence of the method. The residual error infinity norms are given by

$$\begin{array}{c}\hfill \mathrm{R}\left[\mathrm{f}\right]=\left|\right|{\mathsf{\Lambda }}_f{[f,\theta ,\varphi ]\left|\right|}_{\infty },\mathrm{R}\left[\theta \right]=\left|\right|{\mathsf{\Lambda }}_{\theta }{[f,\theta ,\varphi ]\left|\right|}_{\infty },\mathrm{R}\left[\varphi \right]=\left|\right|{\mathsf{\Lambda }}_{\varphi }[f,\theta ,\varphi ]{\left|\right|}_{\infty }.\end{array}$$

(29)

Figure 3 depicts the variation of the residual error infinity norms with the number of iterations. The method achieves an error of about ${10}^{-12}$ after 5 iterations. This shows that the BSQLM is a very accurate method for solving non-linear differential problems.

The effects of the ratio of relaxation to retardation times parameter $\lambda$ on the fluid velocity, temperature and concentration are depicted in Figure 4, Figure 5 and Figure 6. It is observed that enhancing $\lambda$ results in an increase in the fluid velocity near the vertical plate but the velocity decreases towards the free stream. Generally, the fluid temperature and concentration are slightly depressed as $\lambda$ is increased. The ratio of relaxation to retardation times is inversely proportional to the retardation time, which is the time the material require to react to deformation. An increase in $\lambda$ signifies weaker retardation time. Figure 7, Figure 8 and Figure 9 depict the influence of the concentration to thermal buoyancy parameter K on the fluid velocity, temperature, and concentration. From Figure 7, it is noted that the momentum boundary thickens near the plate and then diminishes as the free-stream is approached. Both the temperature and solute concentration distributions are reduced with increasing K as displayed in Figure 8 and Figure 9.

The variation of the fluid velocity, temperature, and concentration with increasing Schmidt number is illustrated in Figure 10, Figure 11 and Figure 12. The dimensionless Schmidt number relates the momentum diffusivity to mass diffusivity. The concentration of the fluid is decreased with the enhanced values of $Sc$. A decrease in the concentration has an effect of suppressing the concentration buoyancy effects which, in turn, will reduce the fluid velocity. However, the opposite trend is realised for the fluid temperature. Figure 13, Figure 14 and Figure 15 reflect the effects of increasing the Deborah number $De$ on the fluid velocity, temperature and concentration profiles, respectively. The importance of the Deborah number in rheology is to define the fluidity of any material under the given specific conditions. A fluid becomes more solid-like with increasing values of the Deborah number, hence a significant decrease in the fluid velocity profile near the vertical plate. For both the fluid temperature and concentration, a similar behaviour is observed when the Deborah number is increased. It is also noted that in the concentration boundary layer, the Deborah number has an insignificant effect on the fluid far away from the vertical plate.

The effects of increasing the Prandtl number on the velocity, temperature, and concentration profiles of the Jeffrey fluid are illustrated in Figure 16, Figure 17 and Figure 18. It is observed that both the fluid velocity and temperature are depressed as $Pr$ increases whilst the concentration is enhanced. The Prandtl number approximates the ratio of the momentum diffusion to thermal diffusion. Increasing the $Pr$ values is equivalent to decreasing thermal diffusion. This will result in the thinning of the thermal boundary layer and hence reduced fluid temperature. Figure 19, Figure 20 and Figure 21 present the influence of the chemical reaction parameter on the fluid velocity, temperature, and concentration profiles, respectively. With increasing $\gamma$, the fluid velocity slightly decreases near the vertical plate and also slightly increase far away from the plate. Increasing the chemical reaction parameter has an effect of reducing the fluid solute concentration. Physically, increasing $\gamma$ results in a reduced concentration boundary layer due to reduced species concentration.

Figure 22, Figure 23 and Figure 24 show that the fluid velocity and temperature profiles are enhanced by an increase in the thermal radiation parameter, $R_d$, whilst the concentration is depressed. Physically, increasing the thermal radiation parameter means more heat is emitted by the Jeffrey fluid causing a rise in the temperature and thickening of the thermal boundary layer. The Eckert number $Ec$ is a dimensionless number that provides a measure of the effects of self-heating of a fluid due to dissipation effects. From Figure 25, Figure 26 and Figure 27, it can be seen that an increasing $Ec$ values enhances velocity and temperature profiles. An opposite trend is observed in the case of the fluid concentration. Physically, increased values of the $Ec$ means more energy dissipating in the thermal boundary layer and this will cause an increase in the fluid temperature.

Figure 28, Figure 29 and Figure 30 depict the variations of the fluid velocity, temperature, and concentration with the magnetic field parameter, M. It is shown that the velocity of the fluid decreases with increasing values of M whilst both the temperature and concentration profiles are enhanced. The velocity and the momentum boundary layer of the fluid are reduced due to the resistive Lorentz force which tends to oppose the fluid flow. The presence of the heat source is seen to cause an increase in the fluid velocity and temperature, as shown in Figure 31 and Figure 32. However, it is revealed in Figure 33 that the concentration profile is suppressed as the values of $\Delta$ are increased. Increasing heat source generates some heat in the thermal boundary layer thickness and the temperature distribution is also increased.

Table 1. Computed values of the skin friction coefficient, local Nusselt number and local Sherwood number for different values of $Pr$, $Sc$, $De$, M and $Ec$ when $\Delta =0.1$, $K=0.1$, $R=0.5$, $\gamma =0.5$, $\xi =1.0$.

$\mathit{Pr}$	$\mathit{Sc}$	$\mathit{\lambda }$	$\mathit{De}$	M	${\mathit{E}}_{\mathit{c}}$	${\mathit{f}}^{″}(0,\mathit{\xi })$	$-{\mathit{\theta }}^{\prime }(0,\mathit{\xi })$	$-{\mathit{\varphi }}^{\prime }(0,\mathit{\xi })$
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
0.2						0.841704	0.434430	2.055455
0.3						0.705153	0.600888	2.036738
0.4						0.601525	0.767872	2.022356
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
	1.0					1.023224	0.266066	3.277721
	1.4					1.019745	0.266026	4.468999
	1.8					1.017726	0.266008	5.659950
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
		0.8				1.333454	0.268194	2.095147
		1.4				1.594027	0.269306	2.104991
		2.0				1.826278	0.270002	2.112276
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
			0.3			0.708783	0.261353	2.058275
			0.5			0.577533	0.258195	2.046650
			0.7			0.500040	0.255768	2.038835
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
				1.0		0.920258	0.259580	2.068256
				1.5		0.828255	0.253782	2.057571
				2.0		0.751381	0.248726	2.048362
0.1	0.6	0.2	0.1	0.5	0.1	1.030290	0.266186	2.080584
					0.5	1.032627	0.260233	2.080870
					0.9	1.034988	0.254216	2.081158
					1.3	1.037374	0.248135	2.081449

displays the variation of the skin friction coefficient, local Nusselt number and the local Sherwood number for various values of the embedded parameters, namely $Pr$, $Sc$, $\lambda$, $De$, M, and $E_c$. It is noted that both the skin friction coefficient and the local Sherwood number are depressed with increasing values of $Pr$, whilst the Nusselt number is increased. When $Pr$ is increased, it means the thermal boundary layer will be thinner relative to the momentum boundary layer. This will result in a large temperature difference and hence an increase in the heat transfer rate. An increase in the Schmidt number causes a decrease in the skin friction coefficient and the Nusselt number whilst the Sherwood number is increased. Increasing the Schmidt number has the effect of decreasing the molecular diffusion and hence a higher concentration of the species. We observe that increasing the values of $\lambda$ increases the skin friction coefficient, local Nusselt number and the local Sherwood number. From Table 1 we can infer that increasing the values of the Deborah number reduces the skin friction coefficient, Nusselt number and Sherwood number. Physically, increasing $De$ means the fluid elasticity is increased whilst the the viscosity is reduced, as well as the frictional force on the plate. We can see that the values of the skin friction coefficient, local Nusselt number, and the local Sherwood number diminish as the values of M are increased. Increasing the Eckert number enhances the skin friction coefficient and the Sherwood number whilst the local Nusselt number is reduced.

5. Conclusions

In this current study, a numerical investigation of the magneto-hydrodynamic laminar boundary flow of a Jeffrey fluid past a vertical plate influenced by viscous dissipation and a heat source or sink was completed. Suitable similarity transformations were used to transform the governing partial differential equations into dimensionless differential equations which were then solved using the bivariate spectral quasi-linearisation method. Some of the most significant findings from this study are the following:

• The bivariate spectral quasi-linearisation method is a very accurate and efficient numerical technique for solving non-linear differential problems;
• Increasing the values of the concentration to thermal buoyancy ratio, radiation parameter, Eckert number, and the heat generation parameter will result in the enhancement of both the fluid velocity and the momentum boundary layer. The opposite trend is observed when the values of the ratio of relaxation to retardation times, Schmidt number, Deborah number, Prandtl number, chemical reaction parameter, and magnetic field parameter are increased;
• The fluid temperature increases with increasing values of the Schmidt number, Deborah number, chemical reaction parameter, radiation parameter, Eckert number, magnetic field parameter, and heat generation parameter, whilst it is reduced with an increase in the values of the ratio of relaxation to retardation times, concentration to thermal buoyancy ratio, and the Prandtl number;
• The fluid concentration is an increasing function of Deborah number, Prandtl number, and the magnetic field parameter whilst decreasing with respect to ratio of relaxation to retardation times, concentration to thermal buoyancy ratio, Schmidt number, chemical reaction parameter, radiation parameter, Eckert number, and the heat generation parameter;
• The local skin friction coefficient is observed to increase with increasing values of the ratio of relaxation to retardation times and the Eckert number. However, it diminishes with increasing values of the Prandtl number, Schmidt number, Deborah number, and magnetic field parameter;
• Increasing the Prandtl number and the ratio of relaxation to retardation times tend to increase the heat transfer rate whilst a decrease is observed when the Schmidt number, Deborah number, magnetic field parameter, and the Eckert number are increased;
• Lastly, the current study observed that the local Sherwood number increase when Schmidt number, the ratio of relaxation to retardation times and the Eckert number are increased whilst it is reduced with increasing Prandtl number, Deborah number, and the magnetic fiels parameter.