Stability Analysis on Magnetohydrodynamic Flow of Casson Fluid over a Shrinking Sheet with Homogeneous-Heterogeneous Reactions

Abstract

Two-dimensional magnetohydrodynamic (MHD) stagnation point flow of incompressible Casson fluid over a shrinking sheet is studied. In the present study, homogeneous-heterogeneous reactions, suction and slip effects are considered. Similarity variables are introduced to transform the governing partial differential equations into non-linear ordinary differential equations. The transformed equations and boundary conditions are then solved using the bvp4c solver in MATLAB. The local skin friction coefficient is tabulated for different values of suction and shrinking parameters. The profiles for fluid velocity and concentration for various parameters are illustrated. It was found that two solutions were obtained at certain ranges of parameters. Then, the bvp4c solver was used to perform stability analysis on the dual solutions. Based on the results, the first solution was more stable and physically meaningful than the other solution. The skin friction coefficient increased when suction increased, but decreased when the magnitude of shrinking parameter increased. Meanwhile, the velocity and concentration profile increased in the presence of a magnetic field. It is also noted that the higher the strength of the homogeneous-heterogeneous reactions, the lower the concentration of reactants.

1. Introduction

Casson fluid is a shear-thinning non-Newtonian fluid that exhibits yield stress, the stress that must be exceeded to make the fluid flow. The fluid acts like a solid and does not flow unless the applied shear stress is larger than the yield stress. Tomato sauce, jelly, honey and even human blood are some of the examples. Many studies have been made on Casson fluid. One of the studies was by Mukhopadhay [1] on the boundary layer flow over a nonlinear stretching surface. Then, this study was extended by Pramanik [2] by considering thermal radiation and an exponentially-stretching surface for the boundary. Recently, the unsteady case was discussed by Khan et al. [3].

The study on the motion of electrically-conducting fluid in the presence of a magnetic field is called magnetohydrodynamics (MHD). In engineering, the application of MHD can be seen in the electromagnetic pump. The pumping motion of this device is caused by the Lorentz force. This force is produced when mutually perpendicular magnetic fields and electric currents are applied in the direction perpendicular to the axis of a pipe containing conducting fluid. Besides that, MHD principles are also used in designing heat exchangers, generators and accelerators. It was found that the strength of magnetic field can affect the fluid flow and temperature distributions inside an enclosure [4]. The study of entropy generation in MHD flow of nanofluid over a rotating disk was done by Rashidi et al. [5]. This study helped to determine the feasibility of using magnetic rotating disk drives in nuclear space propulsion engines. The nanofluid was found to be useful in enhancing the heat transfer rate. Further discussions on the rheological behavior and heat transfer enhancement of nanofluids can be found in the studies by Karimipour et al. [6], Esfe et al. [7], Karimipour et al. [8], Bahrami et al. [9] and Zadkhast et al. [10]. Then, Rashidi et al. [11] extended the study for fluid with variable properties and slip effects. There were also studies on MHD flow of Casson fluid. Nadeem et al. [12] studied the MHD flow with an exponentially-shrinking sheet as the boundary. Raju et al. [13] then discussed the flow over an exponentially-stretching sheet. Thermal radiation, viscous dissipation and chemical reaction were considered in the study. Later, Reddy [14] extended this study by excluding the viscous dissipation effect, and the surface was considered as an inclined stretching surface. The MHD stagnation point flow was discussed by Medikare et al. [15] where the viscous dissipation effect was taken into account in the study. On the other hand, Qing et al. [16] did a study on MHD flow of Casson nanofluid by considering thermal radiation and chemical reactions.

Homogeneous and heterogeneous reactions take place in many chemical reacting systems. Homogeneous reaction occurs at the same phase in bulk of the fluid, whereas heterogeneous reaction occurs at two or more phases commonly on the catalytic surfaces. One of the earliest study was done by Merkin [17]. He studied the isothermal homogeneous-heterogeneous reactions in boundary layer flow over a flat surface. It was found that homogeneous reaction dominates downstream. This leads to other studies related to homogeneous-heterogeneous reactions including the study by Rana et al. [18] on the mixed convection oblique flow of Casson fluid in the presence of internal heating. Then, the study of the reactions on stagnation point flow was done by Raju et al. [19] where the induced magnetic field effect was considered. It was found that the induced magnetic field helped to improve the heat transfer rate. Sheikh and Abbas [20] discussed the reactions on flow over a stretching and shrinking surfaces in the presence of uniform suction and slip effects. Next, Khan et al. [21] expanded the study by Raju et al. [19] into MHD flow with the homogeneous heat effect and omitting the induced magnetic field. The magnetic parameter was found to reduce the velocity. Recently, Khan et al. [22] did a study on a flow with viscous dissipation and Joule heating.

Nowadays, there has been growing interest in the stability analysis of dual solutions. The analysis helps determine the stability and significance of the solutions to the problem. Merkin [23] performed stability analysis on both of the solutions obtained in the study of mixed convection. By considering the unsteady problem, it was found that the upper branch of the solutions was stable while the lower branch was not stable. The same conclusion was obtained in the study by Weidman et al. [24]. Mahapatra and Nandy [25] performed stability analysis on the solutions obtained in stagnation point flow of viscous fluid over a nonlinear shrinking surface. Then, in the study by Roşca and Pop [26], the stability analysis was performed using the bvp4c solver in MATLAB. Next, the analysis was also done by Awaludin et al. [27] on the dual solutions obtained in the shrinking sheet case of stagnation point flow.

Based on the previous research, this study will extend the problem proposed by Sheikh and Abbas [20] to the case of magnetohydrodynamic flow over a shrinking sheet. The differential equations and the boundary conditions will be solved using the bvp4c solver in MATLAB. Stability analysis will then be performed on the dual solutions obtained in this study.

2. Materials and Methods

2.1. Mathematical Formulation

Consider a steady flow of Casson fluid towards a shrinking sheet. Slip effects are considered at the boundary. The Cartesian plane is used to represent the problem such that the shrinking sheet is placed along the x-axis, as shown in Figure 1. The fluid flow is assumed to be confined in the region $y>0$ towards the sheet with $x=0$ as the stagnation point. The sheet is shrunk linearly with a velocity of $u_w\left(x\right)=mx$ where $m<0$ is for the shrinking case and the velocity of external flow is $u_e\left(x\right)=cx$ where $c>0$ is the strength of stagnation flow. A magnetic field with strength $B_0$ is applied in the direction perpendicular to the sheet (y-direction). The induced magnetic field is neglected in this study due to the assumption of a small magnetic Reynolds number.

In the present study, homogeneous-heterogeneous reactions involving two chemical species A and B in the boundary layer flow of Casson fluid near a stagnation point are considered. The problem is analyzed based on a simple model proposed by Merkin [17]. The homogeneous and heterogeneous reactions can be represented by the following equations, respectively,

$$\begin{array}{cc}\hfill A+2B& \to 3B,rate=k_{c}a{b}^2,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill A& \to B,rate=k_{s}a,\hfill \end{array}$$

(2)

where $k_c$ and $k_s$ are the rate constants, while a and b are the respected concentrations of chemical species A and B. The two reaction processes are assumed to be isothermal with a concentration, $a_0$ for reactant A and no autocatalyst B in the ambient fluid.

Based on the stated assumptions, the governing equations are:

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

(3)

$$\begin{array}{cc}\hfill u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}& =u_e\frac{d{u}_e}{dx}+\nu \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\sigma B_0^2}{\rho }(u_e-u),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill u\frac{\partial a}{\partial x}+v\frac{\partial a}{\partial y}& =D_A\frac{{\partial }^{2}a}{\partial y^2}-k_{c}a{b}^2,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill u\frac{\partial b}{\partial x}+v\frac{\partial b}{\partial y}& =D_B\frac{{\partial }^{2}b}{\partial y^2}+k_{c}a{b}^2,\hfill \end{array}$$

(6)

with the boundary conditions:

$$\begin{array}{c}\hfill u\left(0\right)=u_w\left(x\right)=mx+L\left(1+\frac{1}{\beta }\right)\frac{\partial u}{\partial y},v\left(0\right)=-v_w,\\ \hfill D_A\frac{\partial a}{\partial y}{|}_{y=0}=k_{s}a\left(0\right),D_B\frac{\partial b}{\partial y}{|}_{y=0}=-k_{s}a\left(0\right),\\ \hfill u(\infty )=u_e\left(x\right)=cx,a(\infty )=a_0,b(\infty )=0,\end{array}$$

(7)

where u and v are the velocity components in the x- and y-directions, respectively, $\nu$ is the kinematic viscosity, $\sigma$ is the electrical conductivity of the fluid, $\rho$ is the density, $D_A$ and $D_B$ are the respective diffusion coefficients for species A and B, m and c are constants with dimension (time)${}^{-1}$, L is the velocity slip length parameter, $v_w$ is the constant mass transfer velocity and $\beta ={\mu }_B\sqrt{2{\pi }_c}/p_y$ is the Casson parameter where ${\mu }_B$ is the plastic dynamic viscosity of Casson fluid, $p_y$ is the yield stress and ${\pi }_c$ is the critical value for the product of the component of the deformation rate with itself.

The following dimensionless variables are introduced:

$$\begin{array}{c}\hfill \eta =\sqrt{\frac{c}{\nu }}y,u=cx{f}^{\prime }\left(\eta \right),v=-\sqrt{c\nu }f\left(\eta \right),g\left(\eta \right)=\frac{a}{a_0},h\left(\eta \right)=\frac{b}{a_0},\end{array}$$

(8)

for similarity transformations of the governing equations. Substituting (8) into Equations (3)–(7) will result in the the following ordinary differential equations and boundary conditions:

$$\begin{array}{c}\hfill \left(1+\frac{1}{\beta }\right)f^{‴}+f{f}^{″}-f^{\prime 2}+M^2(1-f^{\prime })+1=0,\end{array}$$

(9)

$$\begin{array}{c}\hfill \frac{1}{Sc}{g}^{″}+f{g}^{\prime }-Kg{h}^2=0,\end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{\delta }{Sc}{h}^{″}+f{h}^{\prime }+Kg{h}^2=0,\end{array}$$

(11)

$$\begin{array}{c}\hfill f\left(0\right)=S,f^{\prime }\left(0\right)=\lambda +\upsilon \left(1+\frac{1}{\beta }\right)f^{″}\left(0\right),f^{\prime }(\infty )=1,\end{array}$$

(12)

$$\begin{array}{c}\hfill g^{\prime }\left(0\right)=K_{s}g\left(0\right),g(\infty )=1,\end{array}$$

(13)

$$\begin{array}{c}\hfill \delta h^{\prime }\left(0\right)=-K_{s}g\left(0\right),h(\infty )=0.\end{array}$$

(14)

In the above equations, the prime indicates the ordinary derivative with respect to $\eta$. $M^2=\frac{\sigma B_0^2}{\rho c}$ is the magnetic parameter; $Sc=\frac{\nu }{D_A}$ is the Schmidt number; $K=\frac{k_c{a}_0^2}{c}$ is the strength of homogeneous reaction; while $K_s=\frac{k_{s}R{e}^{-1/2}}{D_A}$ is the strength of heterogeneous reaction with the Reynolds number, $Re=\frac{c}{\nu }$; $\delta =\frac{D_B}{D_A}$ is the ratio of diffusion coefficients; $S=\frac{v_w}{\sqrt{c\nu }}$ is the mass suction parameter with $S>0$; $\lambda =\frac{m}{c}$ is the ratio of the shrinking rate to the rate of external flow with $\lambda <0$ for a shrinking sheet; and $\upsilon =L\sqrt{\frac{c}{\nu }}$ is the velocity slip parameter.

Chaudhary and Merkin [28] stated that the diffusion coefficients of species A and B are of a comparable size in most applications. Thus, it can be assumed that $D_A$ is equal to $D_B$ and $\delta =1$. Therefore, Equations (10) and (11) can be reduced to:

$$\begin{array}{c}\hfill \frac{1}{Sc}{g}^{″}+f{g}^{\prime }-Kg{(1-g)}^2=0,\end{array}$$

(15)

and the boundary conditions are reduced to:

$$\begin{array}{c}\hfill g^{\prime }\left(0\right)=K_{s}g\left(0\right),g(\infty )=1.\end{array}$$

(16)

The transformed equations and boundary conditions are then solved using the bvp4c solver in MATLAB. A detailed derivation of Equations (9)–(16) is provided in Appendix A.

The dimensionless skin friction coefficient $C_f$ is given by:

$$\begin{array}{c}\hfill R{e}_x^{1/2}{C}_f=\left(1+\frac{1}{\beta }\right)f^{″}\left(0\right),\end{array}$$

with $R{e}_x=x{u}_e/\nu$ as the local Reynolds number.

2.2. Stability Analysis

Stability analysis is performed by considering the problem as an unsteady problem. Equations (4)–(6) are replaced by:

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}& =u_e\frac{\partial u_e}{\partial x}+\nu \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\sigma B_0^2}{\rho }(u_e-u),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \frac{\partial a}{\partial t}+u\frac{\partial a}{\partial x}+v\frac{\partial a}{\partial y}& =D_A\frac{{\partial }^{2}a}{\partial y^2}-k_{c}a{b}^2,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \frac{\partial b}{\partial t}+u\frac{\partial b}{\partial x}+v\frac{\partial b}{\partial y}& =D_B\frac{{\partial }^{2}b}{\partial y^2}+k_{c}a{b}^2,\hfill \end{array}$$

(19)

in which t represents time. New dimensionless variables are introduced to replace the variables in (8),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \eta =\sqrt{\frac{c}{\nu }}y,u=cx{f}^{\prime }(\eta ,\tau ),v& =-\sqrt{c\nu }f(\eta ,\tau ),g(\eta ,\tau )=\frac{a}{a_0},\hfill \\ \hfill h(\eta ,\tau )& =\frac{b}{a_0},\tau =ct,\hfill \end{array}\end{array}$$

(20)

so that Equations (9)–(11) can be replaced by:

$$\begin{array}{cc}\hfill \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{3}f}{\partial {\eta }^3}+f\frac{{\partial }^{2}f}{\partial {\eta }^2}-{\left(\frac{\partial f}{\partial \eta }\right)}^2+M^2\left(1-\frac{\partial f}{\partial \eta }\right)+1-\frac{{\partial }^{2}f}{\partial \eta \partial \tau }& =0,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \frac{1}{Sc}\frac{{\partial }^{2}g}{\partial {\eta }^2}+f\frac{\partial g}{\partial \eta }-Kg{h}^2-\frac{\partial g}{\partial \tau }& =0,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{\delta }{Sc}\frac{{\partial }^{2}h}{\partial {\eta }^2}+f\frac{\partial h}{\partial \eta }+Kg{h}^2-\frac{\partial h}{\partial \tau }& =0,\hfill \end{array}$$

(23)

and the boundary conditions become:

$$\begin{array}{cc}\hfill \frac{\partial f}{\partial \eta }(0,\tau )& =\lambda +\upsilon \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{2}f}{\partial {\eta }^2}(0,\tau ),f(0,\tau )=S,\frac{\partial f}{\partial \eta }(\infty ,\tau )=1,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \frac{\partial g}{\partial \eta }(0,\tau )& =K_{s}g(0,\tau ),g(\infty ,\tau )=1,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \delta \frac{\partial h}{\partial \eta }(0,\tau )& =-K_{s}g(0,\tau ),h(\infty ,\tau )=0.\hfill \end{array}$$

(26)

Then, Equations (22) and (23) are reduced to the following equation:

$$\begin{array}{c}\hfill \frac{1}{Sc}\frac{{\partial }^{2}g}{\partial {\eta }^2}+f\frac{\partial g}{\partial \eta }-Kg{(1-g)}^2-\frac{\partial g}{\partial \tau }=0,\end{array}$$

(27)

using the previous assumption where $D_A$ is equal to $D_B$ and $\delta =1$. The boundary conditions are given by:

$$\begin{array}{c}\hfill \frac{\partial g}{\partial \eta }(0,\tau )=K_{s}g(0,\tau ),g(\infty ,\tau )=1.\end{array}$$

(28)

According to Harris et al. [29], in order to test the stability of the steady flow solution $f\left(\eta \right)=f_0\left(\eta \right)$ and $g\left(\eta \right)=g_0\left(\eta \right)$ that satisfies the boundary value problem (9), (12), (15) and (16), we write:

$$\begin{array}{c}\hfill f(\eta ,\tau )=f_0\left(\eta \right)+e^{-\gamma \tau }F(\eta ,\tau ),\\ \hfill g(\eta ,\tau )=g_0\left(\eta \right)+e^{-\gamma \tau }G(\eta ,\tau ),\end{array}$$

(29)

where $F(\eta ,\tau )$ and $G(\eta ,\tau )$ are small relative to $f_0\left(\eta \right)$ and $g_0\left(\eta \right)$ and $\gamma$ is the unknown eigenvalue parameter. Awaludin et al. [27] stated that Equations (21) and (27) along boundary conditions (24) and (28) will yield infinite sets of eigenvalues ${\gamma }_1<{\gamma }_2<{\gamma }_3<{\gamma }_4<\dots$. If the smallest eigenvalue $\gamma <0$, there is initial growth of disturbance in the flow and the flow becomes unstable, whereas if $\gamma >0$, there is initial decay that does not interrupt the flow and the flow is said to be stable. Substituting (29) into Equations (21) and (27) yields the following linearized equations:

$$\begin{array}{cc}\hfill \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{3}F}{\partial {\eta }^3}+f_0\frac{{\partial }^{2}F}{\partial {\eta }^2}+F\frac{{\partial }^2{f}_0}{\partial {\eta }^2}-\left(2\frac{\partial f_0}{\partial \eta }+M^2-\gamma \right)\frac{\partial F}{\partial \eta }-\frac{{\partial }^{2}F}{\partial \eta \partial \tau }& =0,\hfill \end{array}$$

(30)

$$\begin{array}{c}\hfill \frac{1}{Sc}\frac{{\partial }^{2}G}{\partial {\eta }^2}+F\frac{\partial g_0}{\partial \eta }+f_0\frac{\partial G}{\partial \eta }-(K-4K{g}_0+3K{g}_0^2-\gamma )G-\frac{\partial G}{\partial \tau }=0,\end{array}$$

(31)

and the boundary conditions:

$$\begin{array}{c}\hfill F(0,\tau )=0,\frac{\partial F}{\partial \eta }(0,\tau )=\upsilon \left(1+\frac{1}{\beta }\right)\frac{{\partial }^{2}F}{\partial {\eta }^2}(0,\tau ),\frac{\partial G}{\partial \eta }(0,\tau )=K_{s}G(0,\tau ),\end{array}$$

(32)

$$\begin{array}{c}\hfill \frac{\partial F}{\partial \eta }(\eta ,\tau )\to 0,G(\eta ,\tau )\to 0\mathrm{as}\eta \to \infty .\end{array}$$

(33)

Based on Weidman et al. [24], the initial growth or decay in (29) can be identified by setting $\tau$ to 0 so that $F=F_0\left(\eta \right)$ and $G=G_0\left(\eta \right)$ in Equations (30)–(33). Then, the following linear eigenvalue problem is obtained:

$$\begin{array}{cc}\hfill \left(1+\frac{1}{\beta }\right)F_0^{‴}+f_0{F}_0^{″}+f_0^{″}{F}_0-\left(2f_0^{\prime }+M^2-\gamma \right)F_0^{\prime }& =0,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \frac{1}{Sc}{G}_0^{″}+g_0^{\prime }{F}_0+f_0{G}_0^{\prime }-(K-4K{g}_0+3K{g}_0^2-\gamma )G_0& =0,\hfill \end{array}$$

(35)

together with the boundary conditions:

$$\begin{array}{c}\hfill F_0\left(0\right)=0,F_0^{\prime }\left(0\right)=\upsilon \left(1+\frac{1}{\beta }\right)F_0^{″}\left(0\right),G_0^{\prime }\left(0\right)=K_s{G}_0\left(0\right),\end{array}$$

(36)

$$\begin{array}{c}\hfill F_0^{\prime }\left(\eta \right)\to 0,G_0\left(\eta \right)\to 0\mathrm{as}\eta \to \infty .\end{array}$$

(37)

Then, $F_0^{\prime }\left(\eta \right)\to 0$ as $\eta \to \infty$ are chosen to relax. This is because Harris et al. [29] stated that the range of possible eigenvalues can be obtained by relaxing one of the boundary conditions, $F_0\left(\eta \right)$ or $G_0\left(\eta \right)$. Thus, Equations (34) and (35) are solved along the new boundary conditions,

$$\begin{array}{c}\hfill F_0\left(0\right)=0,F_0^{\prime }\left(0\right)=\upsilon \left(1+\frac{1}{\beta }\right)F_0^{″}\left(0\right),G_0^{\prime }\left(0\right)=K_s{G}_0\left(0\right),F_0^{″}\left(0\right)=1,\end{array}$$

(38)

$$\begin{array}{c}\hfill G_0\left(\eta \right)\to 0\mathrm{as}\eta \to \infty .\end{array}$$

(39)

The eigenvalues are computed using the bvp4c solver in MATLAB. The smallest eigenvalue is the eigenvalue with the smallest error. A detailed derivation of Equations (21)–(39) is presented in Appendix B.

2.3. Numerical Method

The nonlinear ordinary differential equations, Equations (9) and (15) subject to the boundary conditions (12) and (16) are solved by using the bvp4c solver in MATLAB. This solver has been widely used by other researchers to solve the boundary value problem. The solver is a finite difference code with fourth order accuracy. In order to use the solver, the equations have to be rewritten as a set of equivalent first order ordinary differential equations. This is done by using the following substitutions:

For Equation (9),

$$\begin{array}{cc}\hfill y\left(1\right)& =f,\hfill \\ \hfill y{\left(1\right)}^{\prime }& =f^{\prime }=y\left(2\right),\hfill \\ \hfill y{\left(2\right)}^{\prime }& =f^{″}=y\left(3\right),\hfill \\ \hfill y{\left(3\right)}^{\prime }& =f^{‴}=\frac{\left[-y\left(1\right)y\left(3\right)+{\left(y\left(2\right)\right)}^2-M^2(1-y\left(2\right))-1\right]}{\left[1+\frac{1}{\beta }\right]}.\hfill \end{array}$$

(40)

For Equation (15),

$$\begin{array}{cc}\hfill y\left(4\right)& =g,\hfill \\ \hfill y{\left(4\right)}^{\prime }& =g^{\prime }=y\left(5\right),\hfill \\ \hfill y{\left(5\right)}^{\prime }& =g^{″}=Sc[-y\left(1\right)y\left(5\right)+Ky\left(4\right){(1-y\left(4\right))}^2].\hfill \end{array}$$

(41)

For the boundary conditions (12) and (16),

$$\begin{array}{cc}\hfill ya\left(1\right)& =S,ya\left(2\right)=\lambda +\upsilon \left(1+\frac{1}{\beta }\right)ya\left(3\right),yb\left(2\right)=1,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill ya\left(5\right)& =K_{s}ya\left(4\right),yb\left(4\right)=1.\hfill \end{array}$$

(43)

Equations (40)–(43) are then coded into the bvp4c solver. The syntax of the solver, sol = bvp4c(@OdeBVP,@OdeBc,solinit,options) contains the function handle @OdeBVP, into which the Equations (40) and (41) are coded. Then, the boundary conditions (42) and (43) are coded into the function handle @OdeBC.

The bvp4c solver applies the collocation formula of three-stage Lobatto IIIa formula. Bilal et al. [30] states that the collocation polynomials provide a $C^1$-continuous solution with fourth order accuracy uniformly in the interval of integration. The interval of integration is set according to the boundary conditions. In this technique, the interval of integration is divided into the subinterval by using a mesh of points. Then, a system of algebraic equations, resulting from the boundary conditions and the collocation conditions imposed on all the subintervals, are solved to obtain the numerical solution. The solver estimates the error of the numerical solution on each of the subintervals. The process is repeated if the solution fails to satisfy the tolerance criteria. The points of the initial mesh and initial approximation of the solution at the mesh points are coded in solinit, while the options is an optional integration argument. The solver is ran, and the results are printed out in the form of numerical solutions and graphs. When the other guess value in the structure solinit results in other solutions that also satisfy the boundary conditions, this indicates that multiple solutions exist to the problem.

The same procedures are done for stability analysis. New substitutions are introduced to rewrite Equations (34) and (35) and the boundary conditions (38) and (39) into first order ordinary differential equations. Let,

$$\begin{array}{ccccc}\hfill y\left(1\right)& =F_0,\hfill & \hfill \mathrm{and}& \hfill s\left(1\right)& =f_0,\hfill \\ \hfill y\left(2\right)& =F_0^{\prime },\hfill & & \hfill s\left(2\right)& =f_0^{\prime },\hfill \\ \hfill y\left(3\right)& =F_0^{″},\hfill & & \hfill s\left(3\right)& =f_0^{″},\hfill \\ \hfill y\left(4\right)& =G_0,\hfill & & \hfill s\left(4\right)& =g_0,\hfill \\ \hfill y\left(5\right)& =G_0^{\prime },\hfill & & \hfill s\left(5\right)& =g_0^{\prime }.\hfill \end{array}$$

The new first order ordinary differential equations are then coded into the solver. The smallest eigenvalue, $\gamma$, is computed by setting an appropriate step size at a command in the solver.

3. Results and Discussion

In the present study, the results on the effects of suction parameter S, shrinking parameter $\lambda$, Casson parameter $\beta$, slip parameter $\upsilon$, magnetic parameter M, the strength of homogeneous reaction parameter K, the strength of heterogeneous reaction parameter $K_s$ and Schmidt number $Sc$ on the fluid velocity and concentration are illustrated in graphs. Besides that, the variation of local skin friction coefficient $f^{″}\left(0\right)$ on various values of S and the effect of S on the concentration gradient $g^{\prime }\left(0\right)$ are also shown in the form of graphs. All the results obtained in this study were of two solutions or also known as the dual solutions. These solutions were obtained at certain ranges of parameters.

Equations (9) and (15) together with the boundary conditions (12) and (16) were solved numerically using the bvp4c solver in MATLAB. The validity of the method used was then checked by comparing the obtained results of $f^{″}\left(0\right)$ with the results from previous studies, as shown in

Table 1. Comparison of the values of local skin friction coefficient, $f^{″}\left(0\right)$, for various values of $\lambda$ when $K=K_s=0.5$, $Sc=1.0$, $M=S=\upsilon =0.0$ and $\beta ={10}^8$.

$\mathit{\lambda }$	Present Study		Kameswaran et al. [31]		Bhattacharyya [32]
	First Solution	Second Solution	First Solution	Second Solution	First Solution	Second Solution
−1.00	1.3288168	0	1.3288169	0	1.3288169	0
−1.15	1.0822311	0.1167021	1.0822312	0.1167021	1.0822316	0.1167023
−1.20	0.9324733	0.2336497	0.9324734	0.2336497	0.9324728	0.2336491
−1.2465	0.5842814	0.5542963	0.5842817	0.5542962	0.5842915	0.5542856

. The computations of $f^{″}\left(0\right)$ in Kameswaran et al. [31] and Bhattacharyya [32] were done by using the bvp4c solver and shooting method, respectively. Thus, the usage of these studies in validating the method used in the present study was suitable. It was found that the results were in good agreement. This reassured that the method used was accurate. In this paper, the computations were done for various values of parameters involved in the flow equations. All the results obtained were in the form of dual solutions. The velocity profiles, concentration profiles, the graph of skin friction and concentration gradient at the surface were plotted.

Figure 2 shows that dual solutions existed when $\lambda >{\lambda }_c$, while no solution was obtained when $\lambda <{\lambda }_c$. ${\lambda }_c$ is the critical point where the solution at this point is unique. The value of $f^{″}\left(0\right)$ for the first solution was observed to increase as S increased, which agrees with the results obtained in

Table 2. Values of the local skin friction coefficient $f^{″}\left(0\right)$ for various values of S and $\lambda$ when $K=0.5$, $K_s=0.5$, $Sc=1.0$, $M=0.2$, $\beta =3.0$ and $\upsilon =0.3$.

S	$\mathit{\lambda }$	${\mathit{f}}^{″}\left(0\right)$
		First Solution	Second Solution
0.0	−1.54	1.098778769	0.680183640
	−1.55	1.066176015	0.725484949
	−1.56	1.020480243	0.783972152
0.1	−1.60	1.201480040	0.720639471
	−1.62	1.140272080	0.808669739
	−1.63	1.090694586	0.871784903
0.3	−1.74	1.401462735	0.841612630
	−1.76	1.352236872	0.920337136
	−1.78	1.269663918	1.032688556
0.5	−1.92	1.576816297	1.060694430
	−1.94	1.515092506	1.154690179
	−1.95	1.465920061	1.220084931

. This shows that the increase in suction caused the wall shear stress to increase. Likewise, the value of $|{\lambda }_c|$ increased as S increased. It is also noted that the value of $f^{″}\left(0\right)$ was always positive. This denotes that the fluid exerted drag force on the sheet.

Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 display the effects of different parameters on the dimensionless fluid velocity, $f^{\prime }\left(\eta \right)$. The values of $f^{\prime }\left(\eta \right)$ are noted to be initially negative. Then, $f^{\prime }\left(\eta \right)$ became positive as $\eta$ increased. The effect of $\lambda$ on $f^{\prime }\left(\eta \right)$ is shown in Figure 3. The increasing value of $\left|\lambda \right|$ caused the value of $f^{\prime }\left(\eta \right)$ to decrease in the first solution and increase in the other solution. The thickness of the momentum boundary layer for the first solution is noted to become larger as $\left|\lambda \right|$ increases. In Figure 4, the increase in $\beta$ caused $f^{\prime }\left(\eta \right)$ to increase and the thickness of momentum boundary layer to decrease in both of the solutions. The large value of $\beta$ indicated that the yield stress was smaller compared to the viscosity and deformation rate. Sheikh and Abbas [20] stated that the yield stress became less when the value of the Casson parameter increased and caused the thickness of the momentum boundary layer to decrease. The presence of tensile stress due to elasticity led to the contraction of boundary layer thickness. Next, Figure 5 shows that the increase in $\upsilon$ caused $f^{\prime }\left(\eta \right)$ of the first solution to increase. The opposite results were obtained for the other solution except for a small part near the surface of the sheet. On the other hand, the effect of S on $f^{\prime }\left(\eta \right)$ is shown in Figure 6. Suction was one of the methods used in delaying the flow separation. Miklavčič and Wang [33] stated that enough suction was needed to maintain the fluid flow on the shrinking surfaces. The presence of suction caused the fluid velocity in the first solution to increase, while the velocity in the other solution reduced. In the first solution, the thickness of the boundary layer reduced as suction increased. The variation of velocity profile with M is illustrated in Figure 7. It can be noted that the higher the value of M, the higher the fluid velocity in the first solution. However, different results were obtained for the second solution where the increase in the value of M caused the velocity to decrease. In the presence of magnetic field, a drag-like force called the Lorentz force is produced. This force provides resistance to the momentum of the fluid particles. This will make the boundary layer thickness to become thinner. In Figure 7, it can be noticed that the thickness of the boundary layer for the first solution becomes thinner as M increased. The boundary layer thickness of the second solution was significantly bigger than the first solution.

The effects of various values of parameters on the concentration, $g\left(\eta \right)$, are illustrated in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15. Figure 8 illustrates the effect of $\lambda$ on $g\left(\eta \right)$. The concentration decreased as $\left|\lambda \right|$ increased in the first solution, whereas the concentration increased as $\left|\lambda \right|$ increased in the second solution. Then, Figure 9 displays that the concentration profiles for both solutions increased as $\beta$ increased. The thickness of concentration boundary layer reduced in both solutions because of the increasing elasticity stress parameter. In Figure 10, the concentration for the first solution increased as slip parameter, $\upsilon$, increased, which was contrary to the second solution. Sheikh and Abbas [20] stated that the concentration increased when the slip parameter increased due to the reduction of available area in the boundary layer, which resulted from the shrinking of the sheet. The same behavior is observed in Figure 11 and Figure 12, where the increase in S and M caused $g\left(\eta \right)$ to increase in the first solution and decrease in the second solution. The thickness of the boundary layer for the first solution became smaller as S and M increased. Besides that, the homogeneous-heterogeneous reactions effects on $g\left(\eta \right)$ are depicted in Figure 13 and Figure 14. It can be noted in both figures that the concentration decreased as K and $K_s$ increased. This was because the reaction rate increased as the values of K and $K_s$ increased, which caused the concentration of the species to decrease. The Schmidt number, $Sc$, is the ratio of the viscous diffusion rate to the molecular diffusion rate. When the value of $Sc>1$, the rate of viscous diffusion was higher than the rate of molecular diffusion. This caused the fluid velocity to decrease and the concentration of reactants to increase, as shown in Figure 15. Last but not least, Figure 16 shows that the presence of S caused the concentration gradient at the surface, $g^{\prime }\left(0\right)$, to increase as $\beta$ increased.

Since the solutions obtained in this paper were of dual solutions, stability analysis was done for the solutions. The smallest eigenvalue, $\gamma$, was calculated by performing the analysis using the bvp4c solver in MATLAB. The positive value of $\gamma$ implied stable flow, while the negative value of $\gamma$ implied unstable flow.

Table 3. Smallest eigenvalues $\gamma$ for different values of $\lambda$ when $K=0.5$, $K_s=0.5$, $Sc=1.0$ and $M=0.2$.

S	$\mathit{\lambda }$	$\mathbf{\gamma }$
		First Solution	Second Solution
0	−1.567	0.11120	−0.10946
	−1.569	0.03000	−0.02988
	−1.5691	0.01824	−0.01819
0.1	−1.635	0.12316	−0.12112
	−1.637	0.05977	−0.05929
	−1.6376	0.01136	−0.01134
0.3	−1.785	0.14267	−0.14016
	−1.788	0.05146	−0.05113
	−1.7884	0.01772	−0.01768
0.5	−1.955	0.14708	−0.14462
	−1.958	0.06012	−0.05971
	−1.9586	0.00710	−0.00710

shows the values of $\gamma$ obtained when $K=K_s=0.5$, $Sc=1.0$ and $M=0.2$. It can be observed that the values of $\gamma$ were approaching zero as $\lambda$ approached the critical value, ${\lambda }_c$. Besides that, it can be noted that $\gamma >0$ for the first solution and $\gamma <0$ for the second solution. Therefore, the first solution is said to be stable and physically meaningful. This agrees with the above discussions. The results of the first solution seemed to be more fitted with the physical explanation provided from previous studies.

4. Conclusions

A study on homogeneous-heterogeneous reactions on MHD flow of Casson fluid was carried out. A shrinking sheet in the presence of suction and slip effects was considered as the boundary. Similarity transformations were introduced to convert the partial differential equations into ordinary differential equations. The bvp4c solver was then used to solve the transformed equations and boundary conditions. The results showed that dual solutions existed at certain ranges of parameters. The results were presented in the form of tables and graphs. The effect of various parameters on fluid velocity and concentration were analyzed and discussed. Since the problem contained two solutions, stability analysis was done to determine the stability of the solutions. The values of the smallest eigenvalue were computed using the bvp4c solver and presented in a table. Based on the results, the first solution was stable and physically meaningful. Then, the following were concluded based on the findings of the first solution:

• The skin friction coefficient increased as suction increased, which indicates that the wall shear stress increased with increasing suction.
• The skin friction coefficient decreased as the magnitude of the shrinking parameter increased.
• The fluid velocity and concentration decreased as the value of the shrinking parameter increased.
• The fluid velocity and concentration increased in the presence of the magnetic parameter.
• The fluid velocity and concentration increased as the value of the Casson parameter, slip parameter and suction increased.
• The stronger the strength of homogeneous-heterogeneous reactions, the lower the concentration.
• The bigger the value of the Schmidt number, the higher the concentration.