Stability analysis of triple solutions of Casson nanofluid past on a vertical exponentially stretching/shrinking sheet

Abstract

The MHD two dimensional boundary layer flow of Casson nanofluid on an exponential stretching/shrinking sheet is considered with effects of radiation parameter, nanoparticles volume fractions (i.e. Fe₃O₄ and Ti₆Al₄V) and thermal convective boundary condition. The partial differential equations are transformed into ordinary differential equations by means of similarity transformations. The solutions of the transferred equations are achieved numerically with the help of shooting technique in Maple software. At different ranges of involved physical parameters, triple solutions are found. Therefore, stability analysis is performed by bvp4c in MATLAB to find the stable and physically reliable solution. Impacts of the physical parameter are presented through graphs and tables. Mainly, it is found that an increase in Casson and suction parameters decrease the corresponding velocity profiles while increase in Prandtl number, stretching/shrinking, and suction parameter decrease the temperature profiles. Furthermore, an increase in nanoparticles volumetric fraction, radiation and magnetic parameters as well as Biot number increase the temperature profiles and their thermal boundary layer thicknesses.

Introduction

The boundary layer flow and heat transfer characteristics of the fluids have wide range of applications in industrial and engineering zones. The processes involving flow on stretching and shrinking sheets (surfaces) have also achieved much importance due to its growing needs in industrial sectors for example, in glass production industries, paper production, polymer extrusion and plastic films, etc. In this regard, first significant attempt in study of stretching surface was reported by Crane. ¹ The Crane’s work was further extended by Chen ² by adding mass transfer characteristics in various physical phenomena. The heat and the mass transfer characteristics on boundary layer fluid flow over exponentially stretching surface was examined by Magyari and Keller. ³ Elbashbeshy ⁴ considered the fluid flow past on the exponentially stretching surface. An investigation on two-dimensional steady boundary layer flow and heat transfer characteristics on exponential stretching surface with thermal radiation effect was considered by Bidin and Nazar. ⁵ Recently, few investigators have examined the problems of the fluid flows over stretching surface (see Khan et al. ⁶ and Zehra et al. ⁷ ).

Whereas, since many years, the heat transfer and thermal conductivity enhancement have become an important field of research. In this regard, most of the researchers have mainly focused on the study of nanofluids which are prepared by suspending nano-sized particles of different solids materials in different conventional fluids that are named as a base fluids of the nanofluids. These fluids can be utilized as a working fluids instead of conventional common fluids in future due to possessing high thermal conductivity, heat capacity, and heat transfer enhancement. In this regard, Choi and Eastman ⁸ were the first who prepared an innovative fluid by suspending solid nano-sized particles in the common base fluid. Lee et al. ⁹ stated that nanofluids possess outstanding heat transfer characteristics as compare to common fluids. Nanofluids can be used in various energetic systems of modern science and technology, like as natural convection enclosures, radiators, cooling of the nuclear systems, etc. Nanofluids are prepared by using Newtonian or non-Newtonian base fluids through suspensions of different types of solid nanoparticles that possess enhanced thermo physical properties as compare to base fluids. To understand the flow and heat transfer characteristics of nanofluids on mathematical viewpoint, Buongiorno ¹⁰ and Tiwari and Das ¹¹ presented mathematical models. Furthermore, the boundary layer nanofluid flow on a nonlinear stretching sheet was studied by Rana and Bhargava. ¹² Mabood et al. ¹³ examined the magnetohydrodynamic (MHD) nanofluid flow and examined the heat transfer characteristics for nonlinear stretching sheet. The studies of the boundary layer nanofluids flows over stretching/shrinking surfaces have been considered by many researchers (see Lund et al. ¹⁴ , Anwar et al. ¹⁵ , Nadeem et al., ¹⁶ and Abbas et al. ¹⁷ ).

In daily life, there are many applications related to fluids such as condensed milk, shampoos, tomatoes juice, printing ink, muds, and paints, etc. All these types of fluids show different types of the characters that are not understandable properly by the Newtonian theory of the fluids. Therefore, it is compulsory to study the non-Newtonians fluids to understand the phenomena of such type of fluids which are full with complexities. The properties related to non-Newtonian fluids may not be stated using a single model of non-Newtonian fluids properly. Therefore, different models have been introduced that are present in literature which are particularly categorized in three different models, namely integral, rate and differential type fluid. In 1959, a model proposed by Casson ¹⁸ that was later named as the Casson fluid flow model in which Casson fluid shows a yield stress. It shows an infinite viscosity at the zero shear stress while it is a type of the shear thinning fluid. In case of the Casson fluid, if applied shear stress is less than yield stress than no flow occurs, whereas, zero viscosity occurs at infinite shear rate, that is, when the shear stress is greater than the applied yield stress, fluid starts moving, whereas fluid show same behavior as a solid if the yield stress is greater than the applied shear stress. The purpose of the present study is to investigate the mixed convection magnetohydrodynamic (MHD) boundary layer flow and the heat transfer characteristics of Fe₃O₄ and Ti₆Al₄V Casson based nanofluid on exponentially stretching/shrinking surface by using Tiwari and Das’s ¹¹ model. Present analysis is also modification of work done by Rohni et al. ¹⁹ and Hafidzuddin et al. ²⁰ The effects of Casson fluid flow are discussed by many researchers (see Hafidzuddin et al. ²⁰ , Amjad et al. ²¹ ).

In this article, similarity transformations are used to transfer partial differential equations into the system of ordinary differential equations that is requirement of the numerical method which is used here. At first, the obtained ordinary differential equations are numerically solved by employing shooting method in Maple software. The obtained results of the equations show that there exist three solutions for each value of the parameters which remains continue at different limits of the pertinent parameters. To check the reliability of the solutions, the stability analysis of the solutions is done by the bvp4c method in MATLAB. The obtained results of the present paper will prove supportive for researchers which are interested to find the impact of different physical fluid parameters on Casson-based nanofluids with the use of the different nanoparticles by applying Tiwari and Das’s model. This study can be beneficial for those which are interested to work on nanofluids on mathematical and experimental basis because in present study two types of the nanoparticles are used and checked their flow and heat transfer characteristics. On other hand, in case of the rising multiple solutions, there has been given process of the stability analysis which indicates that among three obtained solutions, only one solution is reliable for the present study.

Formulation of problem

The two dimensional mixed convection magneto hydro-dynamic (MHD) boundary layer flow and heat transfer of Fe₃O₄ and Ti₆Al₄V-Casson nanofluid over an exponential stretching/shrinking sheet with thermal radiations effect is considered. The Tiwari and Das’s model is used to develop the equations of the present problem. A continuous magnetic field $B(x)$ is used vertically across in the flow direction. The flow geometry and its coordinate system are denoted in the Figure 1. According to Amjad et al., ²² and Ali et al., ²³ the Casson fluid flow rheological equations are considered as:

$${\tau }_{ij}=2\left({\mu }_B+\left(\frac{P_y}{\sqrt{2\pi }}\right)\right)e_{ij},\mathrm{when}\pi >{\pi }_c,$$

(1)

$${\tau }_{ij}=2\left({\mu }_B+\left(\frac{P_y}{\sqrt{2{\pi }_c}}\right)\right)e_{ij},\mathrm{when}\pi <{\pi }_c,$$

(2)

Figure 1.

Physical flow geometry with the coordinate system.

where, ${\tau }_{\mathrm{ij}}$ is the stress tension component, $\pi$ is product of the deformation rate by itself, $e_{\mathrm{ij}}$ is $(i,j)$ th component of deformation rate, ${\pi }_c$ denotes critical with value base on non-Newtonian model and the $P_y(=\frac{{\mu }_B\sqrt{2\pi }}{\beta })$ is the symbol of the yield stress of fluid, and ${\mu }_B$ is a plastic dynamical viscosity of the non-Newtonian fluid. In case of the Casson fluid flow, the dynamic viscosity $(\mu )$ and kinematic viscosity $(\vartheta )$ are defined as,

$$\mu ={\mu }_B+\left(\frac{P_y}{\sqrt{2\pi }}\right)$$

$$\mu ={\mu }_B+\left(\frac{\frac{{\mu }_B\sqrt{2\pi }}{\beta }}{\sqrt{2\pi }}\right)$$

$$\vartheta =\frac{{\mu }_B}{\rho }(1+\frac{1}{\beta })$$

Under above assumptions, the continuity, momentum and the energy equations that defines the Casson nanofluid in accordance with Tiwari and Das model and following to Rohni et al. ¹⁹ can be written as:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$$

(3)

$$\begin{array}{c}u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{1}{{\rho }_{\mathrm{nf}}}\hfill \\ [{\mu }_{\mathrm{nf}}(1+\frac{1}{\beta })\frac{{\partial }^{2}u}{\partial y^2}-\sigma {B_0}^2(x)u+g{\left(\rho {\beta }^{*}\right)}_{\mathrm{nf}}(T-T_{\infty })],\hfill \end{array}$$

(4)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\alpha }_{\mathrm{nf}}\frac{{\partial }^{2}T}{\partial y^2}-\frac{1}{{\left(\rho c_p\right)}_{\mathrm{nf}}}\frac{\partial q_r}{\partial y},$$

(5)

here, u and v represent the components of velocities along the x and y axes. The surface is assumed vertically stretched with velocity $u_w=U_w{e}^{\frac{x}{L}}$ , here $U_w$ is constant of velocity and $L$ is reference length. The surface temperature is taken as $T_w=T_{\infty }+T_0{e}^{\frac{2x}{L}}$ where $T_0$ is constants and $T_{\infty }$ is ambient temperature. ${\rho }_{\mathrm{nf}}=(1-\varphi ){\rho }_f+\varphi {\rho }_s$ is the effective density, ${\left(\rho {\beta }^{*}\right)}_{\mathrm{nf}}=(1-\varphi ){\left(\rho {\beta }^{*}\right)}_f+\varphi {\left(\rho {\beta }^{*}\right)}_s$ is effective thermal expansion coefficient of nanofluid, ${\mu }_{\mathrm{nf}}=\frac{{\mu }_f}{{(1-\varphi )}^{2.5}}$ is effective dynamic viscosity, ${\left(\rho c_p\right)}_{\mathrm{nf}}=(1-\varphi ){\left(\rho c_p\right)}_f+\varphi {\left(\rho c_p\right)}_s$ is effective heat capacitance, ${\alpha }_{\mathrm{nf}}=\frac{k_{\mathrm{nf}}}{{\left(\rho c_p\right)}_{\mathrm{nf}}}$ is nanofluid thermal diffusivity and $k_{\mathrm{nf}}=\frac{(k_s+2k_f)-2\varphi (k_f-k_s)}{(k_s+2k_f)+\varphi (k_f-k_s)}{k}_f$ is the nanofluid thermal conductivity. Furthermore, $\lambda$ stands for stretching and shrinking parameter. While by Roseland approximation for the radiation, the radiative heat flux is written as:

$$q_r=\frac{4{\sigma }^{*}}{3k^{*}}\frac{\partial T^4}{\partial y},$$

(6)

here, $k^{*}$ denotes coefficient of mean absorption while Stefan Boltzmann constant is denoted by symbol ${\sigma }^{*}$ . The difference in temperature in flows are supposed to be sufficiently smaller so as $T^4$ should be stated as the linear function of $T$ , applying the truncated form of Taylor series on ambient temperature $T_{\infty }$ and then ignoring the higher order terms due to smaller values, we get:

$$T^4\cong 4T_{\infty }^{3}T-3T^4$$

(7)

The associated boundary conditions are,

$$\begin{array}{c}v=v_w;u=\lambda u_w\left(x\right);-k_f\frac{\partial T}{\partial y}=h_f(T_f-T);\mathrm{at}y=0\\ u\to 0;T\to T_{\infty };C\to C_{\infty }\mathrm{as}y\to \infty .\end{array}$$

(8)

The required similarity solutions are obtained with the help of similarity transformations which are

$$\begin{array}{c}\multicolumn{1}{c}{\psi =\sqrt{2\vartheta {\mathrm{lU}}_w}{e}^{x/2l}f\left(\eta \right),\theta \left(\eta \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }},}\\ \eta =y\sqrt{\frac{U_w}{2\vartheta l}}{e}^{x/2l},\hfill \end{array}$$

(9)

In form of the velocity components, the stream function $\psi$ is

$$u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x},$$

and magnetic field $B(x)$ with constant $B_0$ for exponential surface is assumed to be

$$B=B_0{e}^{x/L},$$

With the help of the transformations given in equation (9), the system of the equations (4), (5), and (8) is written as

$$\begin{array}{l}\left(1+\frac{1}{\beta }\right)f^{\prime \prime \prime }+{\left(1-\varphi \right)}^{2.5}\\ \left[\left(\left(1-\varphi \right)+\varphi \left(\raisebox{1ex}{${\rho }_s$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.\right)\right)\left(f{f}^{\prime \prime }-2f^{\prime 2}\right)-2Mf\prime +\xi \left(\left(1-\varphi \right)+\varphi \frac{{\left(\rho {\beta }^{*}\right)}_s}{{\left(\rho {\beta }^{*}\right)}_f}\right)\theta \right]=0,\end{array}$$

(10)

$$\begin{array}{c}\frac{k_{\mathrm{nf}}/k_f}{\mathrm{Pr}}(1+\frac{4}{3}\mathrm{Rd}){\theta }^{″}+((1-\varphi )+\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})\hfill \\ (f{\theta }^{\prime }-4f^{\prime }\theta )=0,\hfill \end{array}$$

(11)

the boundary conditions are

$$f\left(0\right)=S;f^{\prime }\left(0\right)=\lambda ;\theta \left(0\right)=1+\frac{1}{\mathrm{Bi}}\theta \prime \left(0\right),\mathrm{at}\eta =0,$$

$$f^{\prime }\left(\eta \right)\to 0;\theta \left(\eta \right)\to 0;\mathrm{as}\eta \to \infty .$$

(12)

Moreover, $\xi =\frac{\mathrm{Gr}}{R{e_x}^2}$ is mixed convection parameter where $\mathrm{Gr}=\frac{L^{3}g{{\beta }^{*}}_f{T}_0}{{{\vartheta }_f}^2}$ is Grashof number, $R{e}_x=\frac{L{U}_w}{{\vartheta }_f}$ , $\mathrm{Rd}=\frac{4T_{\infty }^3{\sigma }^{*}}{k_f{k^{*}}_f}$ is radiation parameter $,\mathrm{Pr}=\frac{{\vartheta }_f{\left(\rho c_p\right)}_f}{k_f}$ is the Prandtl number, $M=\frac{L\sigma B^2(x)}{u_w}$ is a magnetic parameter and $\mathrm{S}=-\sqrt{\frac{2L}{\vartheta U_w}}$ is suction parameter, $\lambda$ is a stretching/shrinking parameter, $\mathrm{Bi}=\frac{h_f\sqrt{\frac{2L{\vartheta }_f}{U_w}}}{k_f{e}^{\frac{x}{2L}}}$ is Biot number. The skin friction coefficient $(C_f)$ and the local Nusselt number $(N{u}_x)$ both are considerable quantities which are take as,

$$C_f=\frac{{\tau }_w}{\rho {\left(u_w\right)}^2},N{u}_x=-\frac{L{q}_w}{k_f(T_w-T_{\infty })},$$

(13)

here, ${\tau }_w$ is the shear stress of wall and $q_w$ is heat flux from surface to fluid that can be expressed as

$${\tau }_w={\mu }_{\mathrm{nf}}(1+\frac{1}{\beta }){\left(\frac{\partial u}{\partial y}\right)}_{y=0}\mathrm{and}{q}_w=k_{\mathrm{nf}}{\left(\frac{\partial T}{\partial y}\right)}_{y=0}.$$

(14)

Using equations (9) and (14) in equations (13), we get,

$$\begin{array}{c}{C}_f{\left(2R{e}_x\right)}^{\frac{1}{2}}=\frac{1}{{(1-\varphi )}^{2.5}}(1+\frac{1}{\beta })f^{″}\left(0\right),\hfill \\ N{u}_x{(2/R{e}_x)}^{\frac{1}{2}}=-\frac{k_{\mathrm{nf}}}{k_f}{\theta }^{\prime }\left(0\right).\hfill \end{array}$$

(15)

where, $R{e}_x=\frac{L{u}_w}{{\vartheta }_f}$ is a Reynolds number.

Stability analysis

There occurs triple solutions in this problem and in such cases, it is necessary to perform the stability analysis of the solutions to check the feasibility of the solutions that which one solution is stable and physically reliable. For stability analysis of the solutions, at first the momentum equation (4) and heat transfer equation (5) are written into the unsteady form as done by Zehra et al. ⁷ by introducing new time dependent variable $\tau .$

$$\begin{array}{c}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{1}{{\rho }_{\mathrm{nf}}}\hfill \\ [{\mu }_{\mathrm{nf}}(1+\frac{1}{\beta })\frac{{\partial }^{2}u}{\partial y^2}-\sigma B^2(x)u+g{\left(\rho {\beta }^{*}\right)}_{\mathrm{nf}}(T-T_{\infty })],\hfill \end{array}$$

(16)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{k_{\mathrm{nf}}}{{\left(\rho C_p\right)}_{\mathrm{nf}}}(1+\frac{4\mathrm{Rd}}{3})\frac{{\partial }^{2}T}{\partial y^2},$$

(17)

here, t is time. $\tau$ is taken as a new dimensionless variable, by variable $\tau$ , the equation (9) can be expressed as:

$$\begin{array}{c}\psi =\sqrt{2U_w\vartheta L}{e}^{\frac{x}{2L}}f(\eta ,\tau );\theta (\eta ,\tau )=\frac{T-T_{\infty }}{T_w-T_{\infty }};\hfill \\ \eta =y\sqrt{\frac{U_w}{2L\vartheta }}{e}^{\frac{x}{2L}}\mathrm{and}\tau =\frac{U_w}{2l}{e}^{x/\mathrm{Lt}}.\hfill \end{array}$$

(18)

where, $u=U_w{e}^{\frac{x}{L}}{f}^{\text{'}}(\eta ,\tau ),v=-\sqrt{\frac{U_w{\vartheta }_f}{2L}}{e}^{\frac{x}{2L}}[f(\eta ,\tau )+\eta f^{\text{'}}(\eta ,\tau )]$ . By applying equation (18) in equations (16) and (17), it is obtained as:

$$\begin{array}{l}\left(1+\frac{1}{\beta }\right)\frac{{\partial }^{3}f\left(\eta ,\tau \right)}{\partial {\eta }^3}+{\left(1-\varphi \right)}^{2.5}\\ \left[\begin{array}{l}\left(\left(1-\varphi \right)+\varphi \left(\raisebox{1ex}{${\rho }_s$}\!\left/ \!\raisebox{-1ex}{${\rho }_f$}\right.\right)\right)\left\{\begin{array}{l}f\left(\eta ,\tau \right)\frac{{\partial }^{2}f\left(\eta ,\tau \right)}{\partial {\eta }^2}-2{\left(\frac{\partial f\left(\eta ,\tau \right)}{\partial \eta }\right)}^2\\ -2\tau \frac{{\partial }^{2}f\left(\eta ,\tau \right)}{\partial \eta \partial \tau }\frac{\partial f\left(\eta ,\tau \right)}{\partial \eta }-\frac{{\partial }^{2}f\left(\eta ,\tau \right)}{\partial \eta \partial \tau }\end{array}\right\}\\ -\frac{2L\sigma B^2(x)}{u_w}\frac{\partial f\left(\eta ,\tau \right)}{\partial \eta }+2\xi \left(1-\varphi +\varphi \frac{{\left(\rho {\beta }^{*}\right)}_s}{{\left(\rho {\beta }^{*}\right)}_f}\right)\theta \left(\eta ,\tau \right)\end{array}\right]=0,\end{array}$$

(19)

$$\begin{array}{c}\frac{1}{\mathrm{Pr}}\frac{k_{\mathrm{nf}}}{k_f}(1+\frac{4\mathrm{Rd}}{3})\frac{{\partial }^2\theta (\eta ,\tau )}{\partial {\eta }^2}+(1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})\hfill \\ \left[\begin{array}{c}f(\eta ,\tau )\frac{\partial \theta (\eta ,\tau )}{\partial \eta }\hfill \\ \multicolumn{1}{c}{-4\frac{\partial f(\eta ,\tau )}{\partial \eta }\theta (\eta ,\tau )-2\tau f^{\prime }(\eta ,\tau )\frac{\partial \theta (\eta ,\tau )}{\partial \tau }-\frac{\partial \theta (\eta ,\tau )}{\partial \tau }}\end{array}\right],\hfill \end{array}$$

(20)

boundary conditions equation (12) can be written as,

$$\begin{array}{c}f(0,\tau )=S,\frac{\partial f(0,\tau )}{\partial \eta }=-1,\frac{\partial \theta (0,\tau )}{\partial \eta }=-\mathrm{Bi}(1-\theta (0,\tau )),\\ \frac{\partial f(\eta ,\tau )}{\partial \eta }\to 0,\theta (\eta ,\tau )\to 0\mathrm{as}\eta \to \infty .\end{array}$$

(21)

The perturbed basic solutions $f\left(\eta \right)=f_0\left(\eta \right)$ and $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$ with disturbances are expressed as,

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

(22)

where, $\gamma$ shows the smallest eigenvalue, $F\left(\eta \right)$ and $G\left(\eta \right)$ both are small relative to $f_0\left(\eta \right)$ and ${\theta }_0\left(\eta \right)$ respectively. In the last, we get a system of the linearized eigenvalue problem,

$$\begin{array}{c}\multicolumn{1}{c}{(1+\frac{1}{\beta })F_0^{″\text{'}}+{(1-\varphi )}^{2.5}\left\{\begin{array}{c}((1-\varphi )+\varphi \left(\frac{{\rho }_s}{{\rho }_f}\right))(f_0{F}_0^{″}+F_0{f}_0^{″}-4f_0^{\text{'}}{F}_0^{\text{'}}+2\tau \gamma f_0^{\text{'}}{F}_0^{\text{'}}+\gamma F_0^{\text{'}})\hfill \\ \multicolumn{1}{c}{-2{\mathrm{MF}}_0^{\text{'}}+\xi (1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})G}\end{array}\right\}=0,}\end{array}$$

(23)

$$\frac{1}{\mathrm{Pr}}\frac{k_{\mathrm{nf}}}{k_f}(1+\frac{4\mathrm{Rd}}{3})G_0^{″}+(1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})\{f_0{G}_0^{\text{'}}+F_0{\theta }_0^{\text{'}}-4{f^{\text{'}}}_0{G}_0-4F\text{'}{\theta }_0+2\tau \gamma f{\text{'}}_0{G}_0+\gamma G_0\}=0,$$

(24)

By setting $\tau =0,$ we get,

$$(1+\frac{1}{\beta })F_0^{″\prime }+{(1-\varphi )}^{2.5}\left\{\begin{array}{c}((1-\varphi )+\varphi \left(\frac{{\rho }_s}{{\rho }_f}\right))(f_0{F}_0^{″}+F_0{f}_0^{″}-4f_0^{\prime }{F}_0^{\prime }+\gamma F_0^{\prime })\hfill \\ \multicolumn{1}{c}{-2{\mathrm{MF}}_0^{\prime }+\xi (1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})G}\end{array}\right\}=0,$$

(25)

$$\begin{array}{c}\frac{1}{\mathrm{Pr}}\frac{k_{\mathrm{nf}}}{k_f}(1+\frac{4\mathrm{Rd}}{3})G_0^{″}+(1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})\hfill \\ \{f_0{G}_0^{\prime }+F_0{\theta }_0^{\prime }-4{f^{\prime }}_0{G}_0-4F\prime {\theta }_0+\gamma G_0\}=0,\hfill \end{array}$$

(26)

with boundary conditions,

$$F_0\left(0\right)=0,F_0^{\prime }\left(0\right)=0,G_0\left(0\right)=0,$$

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

(27)

The above linearized equations (25) and (26) with boundary conditions equation ( $27$ ) are need to be solved. Furthermore, to find the smallest eigenvalues ( $\gamma$ ) by using the bvp4c solver in Matlab software. To obtain the smallest eigenvalues, named as $\gamma$ , one of the boundary condition is transferred in the from of the initial condition that was proposed by Harris et al. ²⁴ Here, in this problem, $F_0\prime (\eta )\to 0$ as $\eta \to \infty$ is converted to $F{\prime }_0^{\prime }(0)=1$ . It is mentioned that the smallest obtained negative values of $\gamma$ show the disturbance in initial growth, so the solution is said to be unstable. While, if obtained smallest value of $\gamma$ is positive, then the solution is said to be stable and physically feasible.

Numerical Method

The boundary value problem (BVPs) expressed in equations (10) and (11) subjecting to initial and the boundary conditions specified in equation (12) are solved by shooting technique which is explained by Meade et al. ²⁵ A short procedure is defined as,

$$f^{\prime }=F_P,f^{″}=F_{\mathrm{PP}},{\theta }^{\prime }={\theta }_P,{\varphi }^{\prime }={\varphi }_P,$$

(28)

$$\begin{array}{c}(1+\frac{1}{\beta })F_{\mathrm{PP}}^{\prime }+{(1-\varphi )}^{2.5}\hfill \\ \left[\begin{array}{c}\left\{(1-\varphi +\varphi \left(\frac{{\rho }_s}{{\rho }_f}\right))(F{F}_{\mathrm{PP}}-2{\left(F_P\right)}^2)\right\}\hfill \\ \multicolumn{1}{c}{-2\mathrm{Mf}\prime +2\xi (1-\varphi +\varphi \frac{{\left(\rho {\beta }^{*}\right)}_s}{{\left(\rho {\beta }^{*}\right)}_f})\theta }\end{array}\right]=0,\hfill \end{array}$$

(29)

$$\begin{array}{c}(1+\frac{4\mathrm{Rd}}{3}){\theta }_P^{\prime }+\mathrm{Pr}\frac{k_f}{k_{\mathrm{nf}}}(1-\varphi +\varphi \frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f})\hfill \\ (F{\theta }_P-4F_P\theta )=0,\hfill \end{array}$$

(30)

along the conditions

$$\begin{array}{c}\multicolumn{1}{c}{F\left(0\right)=S,F_P\left(0\right)=\lambda ,\theta \left(0\right)=1+\frac{1}{\mathrm{Bi}}{\theta }_P\left(0\right),}\\ \multicolumn{1}{c}{F_{\mathrm{PP}}\left(0\right)={\alpha }_1,{\theta }_P\left(0\right)={\alpha }_2,}\end{array}$$

(31)

where, ${\alpha }_1$ and ${\alpha }_2$ are unknown initial conditions. Hence, the initial values of ${\alpha }_1$ , ${\alpha }_2$ is crucial that the solution might satisfy the boundary conditions $F_P\left(\eta \right)\to 0,\theta \left(\eta \right)\to 0\mathrm{as}\eta \to \infty$ of the present boundary value problem. The calculation is done by numerical technique named as shooting method in the Maple software. In present study, the boundary layer thickness $({\eta }_{\infty })$ is taken between $5$ and $6$ which provides a good convergence in solutions. The ranges of the velocity boundary layers thicknesses remained comparatively higher than thermal boundary layer thicknesses. Consequently, that provides three different values for $f^{\prime }\prime \left(0\right)$ and $-\theta \prime \left(0\right),$ respectively.

Results and discussion

The equations (10) and (11) that subject to the boundary conditions in equation (12) are numerically solved to investigate the influence of different applied parameters on magnetohydrodynamic (MHD) mixed convection boundary layer Fe₃O₄ and Ti₆Al₄V-Casson nanofluid flow on a vertical exponentially stretching/shrinking sheet. The numerical results of these equations are calculated by means of shooting technique in the Maple software. The values of skin friction coefficient $(f^{\prime }\prime (0)),$ local Nusselt number $(-\theta \prime (0))$ , velocity profiles $(f^{\prime }\left(\eta \right))$ , and the temperature profiles $(\theta (\eta ))$ are examined by taking various values of the physical parameters. To check the validation of our obtained results, the comparison has been made with previously obtained results by Hafidzuddin et al. ²⁰ which is presented in Table 1. The values of the thermo-physical properties of Fe₃O₄ and Ti₆Al₄V nanoparticles as well as Casson fluid are presented in Table 2. Due to occurrence of the triple solutions, there has been done stability analysis to find the feasibility of the solutions that to be checked which one solution is stable and physically feasible. To obtain the stability analysis, we have used another technique that is three stage Lobatto III a formula. This formula is constructed in bvp4c by utilizing the finite difference code. Later on, the bvp4c solver functions are applied to find the results. The obtained smallest eigenvalues are presented in Table 3, which indicates the stability of first and physically feasible solutions due to having positive smallest eigenvalues. While, second and third solutions are found unstable that can not be physically feasible solutions due to having negative smallest eigenvalues.

Table 1.

Comparative of the values of $f″(0)$ and $-\theta \text{'}(0)$ with those of Hafidzuddin et al. ²⁰ for different $S$ when $\varphi =\mathrm{Rd}=\xi =0$ , $\lambda =1$ , $\beta \to \infty ,\mathrm{Bi}\to \infty$ , and $\mathrm{Pr}=0.72$ .

$S$	Hafidzuddin et al. 19		$\mathrm{Present}$
	$f″(0)$	$-\theta \text{'}(0)$	$f″(0)$	$-\theta \text{'}(0)$
0	1.28182	0.767669	1.28179	0.767769
0.6	1.59824	1.014570	1.59825	1.014571

Table 2.

The thermo physical characteristics of Casson fluid and nanoparticles (Hatami and Ganji, ²⁶ Zaib et al. ²⁷ ).

	$\rho (\mathrm{kg}/m^3)$	$C_p(J/\mathrm{Kgk})$	$k(W/\mathrm{mk})$	$\beta \times {10}^{-5}(1/k)$
Casson fluid (C6H9NaO7)	989	4175	0.6376	23
Ti6Al4V	4420	0.56	7.2	5.8
Fe3O4	5180	670	80.4	20.6

Table 3.

The smallest eigenvalues of $\gamma$ in relation with the parameters $M$ and $\beta$ .

$\beta$	$M$	$\gamma$
		First solution	Second solution	Third solution
2	0.5	3.5105	−3.5401	−3.3510
2.3	0.5	4.9318	−3.4763	−2.8452
2	0.3	4.8563	−3.5680	−2.5437
2	0.6	4.8793	−3.4824	−3.5234

Figures 2 to 5 show the variation in the skin friction coefficients $(f″\left(0\right))$ and the local Nusselt number $(-\theta \prime (0))$ along the stretching/shrinking parameter $(\lambda )$ at different values of suction parameter $(S),$ where values of other parameters are kept constant that are shown in each graph of Fe₃O₄ and Ti₆Al₄V Casson nanofluid, respectively. These figures show that triple solutions exist for the both stretching (λ > 0) and shrinking cases ( $\lambda <0$ ). This study show that two solution merge into one another at a critical points that are mentioned as ${\lambda }_{c_0},{\lambda }_{c_1},{\lambda }_{c_2}$ for $S=3,3.5$ , and $4$ , respectively. Furthermore, Figures 2 to 5 illustrate the influence of suction parameter $\left(S\right)$ on the coefficient of skin friction $f^{\prime }\prime \left(0\right)$ and the local Nusselt number $(-\theta \prime \left(0\right))$ of Fe₃O₄ and Ti₆Al₄V-Casson nanofluid due to change in stretching/shrinking parameter $(\lambda )$ , while values of remaining parameters are fixed. The triple solutions are found for $\lambda >{\lambda }_c$ ( ${\lambda }_{c_0},{\lambda }_{c_1}\mathrm{and}{\lambda }_{c_2}$ ), where second and the third solutions join to gather at the critical points ${\lambda }_c$ and first solution remains continuous to exist for stretching and the shrinking cases of $\lambda$ . Moreover, Figures 2 and 3 indicates the effects of $S$ on the skin friction coefficients $(f^{\prime }\prime (0))$ of Fe₃O₄ and Ti₆Al₄V-Casson nanofluid with variation of the $\lambda$ , respectively. These figures denote that the magnitude of $f^{\prime }\prime \left(0\right)$ reduces in the case of stretching surface $(\lambda >0)$ while increases for shrinking case $(\lambda <0)$ in first solution of Fe₃O₄ and Ti₆Al₄V-Casson nanofluid at increasing rate of suction. An increase in $S$ reduces the rate of skin friction $(f^{\prime }\prime (0))$ in second and third solutions. Figures 4 and 5 show the effect of $S$ on Nusselt number $(-\theta (0))$ with variation of $\lambda$ for Fe₃O₄ and Ti₆Al₄V-Casson nanofluids. These figures show that the increasing rate of the $S$ causes to rise in the heat transfer rate $(-\theta \prime (0))$ for first solutions throughout the flow. While the magnitude of $-\theta \prime \left(0\right)$ increases in second and third solutions when $S$ is increased along the variation of $\lambda$ on Fe₃O₄ and Ti₆Al₄V-Casson nanofluid, respectively. Figures 6 and 7 illustrate the comparative result of skin friction coefficient and the Nusselt number related to variation of volume fraction $(\varphi )$ for Fe₃O₄ and Ti₆Al₄V nanoparticles on Casson nanofluid flow. Figure 6 illustrates the suspension of Ti₆Al₄V nanoparticles in Casson fluid provides more drag force as compare to Fe₃O₄ nanoparticles. While Fe₃O₄ nanoparticles provide more rate of heat transfer as compare to Ti₆Al₄V nanoparticles in Casson nanofluid flow that is illustrated in Figure 7. Figure 8 shows the temperature profile of Ti₆Al₄V-Casson nanofluid flow is slightly remains greater than Fe₃O₄-Casson nanofluid flow in first and second solutions. While in the third solution converse result is observed. Furthermore, to avoid the repetition, the velocity and temperature profiles of Fe₃O₄ Casson nanofluid flow are described only due to possessing same trend of results of different physical parameters on the Ti₆Al₄V-Casson nanofluid flow. Where, Figures 9 and 10 display the impact of the volume fractions of Fe₃O₄ nanoparticles on velocity profile and the temperature profile of the Casson nanofluid. Both these profiles are increasing by increasing nanoparticle volume fraction in Casson based nanofluid in first solution clearly, while in the second and the third solutions at start both profiles are increasing but after a point both are decreasing, respectively. Figure 11 shows the variation in velocity profile of Fe₃O₄-Casson nanofluid by increase in the Casson parameter $(\beta )$ . In the first (stable) solution, an increase in $\beta$ decreases velocity profile with its boundary layer thickness throughout the flow. While in second and the third solutions, velocity profiles are decreasing initially, but after a point they are increasing in both solutions. Actually, an increase in the $\beta$ rises the dynamic viscosity. The yield stress that drags the nanofluid to the stretching surface which develops the resistance in flow due to rising the dynamic viscosity. In result, the velocity and the their boundary layer thicknesses are decreasing when parameter $\beta$ is increased. Figures 7 to 12 show the effect of Casson parameter $(\beta )$ on the temperature profile of Fe₃O₄-Casson based nanofluid. From this Figure it can be seemed temperature and the thickness of the thermal boundary layers are increased with increase in $(\beta )$ in first solution. While second and third solution show decreasing behavior of temperature. Furthermore, Figures 13 and 14 reveal that the stretching/shrinking parameter $(\lambda )$ variates the velocity and the temperature profile of Fe₃O₄-Casson nanofluid flow, respectively. Figure 13 demonstrates that an increase in $\lambda ,$ the velocity, and the boundary layers thickness are increased in first (stable) solution clearly. While in the second and in third solutions the fluctuations are observed where increasing and decreasing behaviors are observed. Figure 14 indicates the temperature profile and thermal boundary layer thickness decrease for larger values of $\lambda$ in first (stable) and third (unstable) solutions throughout the flow, while in second solution, it increases. The change in velocity and temperature profiles of Fe₃O₄-Casson nanofluid flow due to the suction parameter $(S)$ is given in Figures 15 and 16, respectively. Figure 15 indicates the velocity profile and its boundary layer thickness decreases with increase in $S$ in first (stable) solution. Actually, flow of the nanofluid comes nearer to solid sheet and velocity boundary layer becomes thinner when $S$ is increased. While in second and third solutions, the velocity profiles decrease initially, but after a point they are increasing in both unstable solutions. It is noted from Figure 16, the temperature profile and the thermal boundary layers decrease when the rate of the suction is increased in all three obtained solutions. This is due to that the flow comes closer to boundary layer by increase in mass suction, thus, thickness of velocity boundary layer takes the thinner form and also temperature losses easily with flow outward. Figure 17 indicates the relation of magnetic parameter $(M)$ with the temperature profile of Fe₃O₄-Casson nanofluid flow. This figure demonstrates the velocity in first and third solutions increases throughout the flow with increase in magnetic parameter and same behaviours can be seen for the thickness of its momentum boundary layer in this particular case. Actually, an increase in magnetic parameter develops a resisting force (Lorentz force) in flow of the fluid, that acts in contrast to the flow direction. This resistance in flow becomes a main reason of enhancement in temperature profile when magnetic parameter is increased. While the second solution indicates the decreasing trend of the temperature profile with its boundary layer thickness. Figure 18 is drawn to express the effect of radiation parameter $(\mathrm{Rd})$ on temperature profile. It can be seen from this figure, the temperature profile and thickness of its thermal boundary layers are increasing in first solution with increase in $\mathrm{Rd}$ . While in the second solution, temperature profile and its boundary layer thickness decrease throughout the flow and in third solution, it increases at the start but after a while it decreases as $\mathrm{Rd}$ is increased. The rise in radiation parameter emit heat energy toward the flow of the fluid, therefore, the temperature field is increased by heat energy. Thus, the high values of the radiation develops thickness of thermal boundary layer. The impact of Biot number $(\mathrm{Bi})$ on temperature profile is illustrated in Figure 19. It is clear that the increasing value of the $\mathrm{Bi}$ increases the temperature profile in first solution and its thermal boundary layer thickness. While in second and third solutions it increases at the beginning but after a while it decreases. Figure 20 is designed to investigate the impact of Prandtl number $\mathrm{Pr}$ on the temperature profile of Fe₃O₄-Casson nanofluid flow. It is concluded that the temperature and its thermal boundary layer thickness are reducing at increasing value of the $\mathrm{Pr}$ in first solution throughout the flow. This is because of the fluids at high value of the $\mathrm{Pr}$ possess the weak thermal diffusivity results decrease of temperature of flowing fluid. Consequently, temperature and thermal boundary layer are decreasing by increasing value of $\mathrm{Pr}$ . Whereas, second and the third solution show fluctuation, initially temperature decreases in both solutions, but after a point, it increases in related solutions.

Figure 2.

Skin friction coefficient of Ti₆Al₄V-Casson nanofluid with variation of $\lambda$ for the different values of $S.$

Figure 3.

Skin friction coefficient of Fe₃O₄-Casson nanofluid with variation of $\lambda$ for the various values of $S.$ .

Figure 4.

Local Nusselt number of Ti₆Al₄V-Casson nanofluid with variation of $\lambda$ for the different values of $S.$

Figure 5.

Local Nusselt number of Fe₃O₄-Casson nanofluid with variation of $\lambda$ for the different values of $S.$ .

Figure 6.

The skin friction coefficient variations with $\varphi$ for the used specific nanoparticles.

Figure 7.

The local Nusselt number variations with $\varphi$ for used specific nanoparticles.

Figure 8.

Temperature profiles for Fe₃O₄ and Ti₆Al₄V Casson based nanofluid.

Figure 9.

Velocity profiles for the various values of $\varphi$ .

Figure 10.

Temperature profiles for the various values of $\varphi$ .

Figure 11.

Velocity profiles for the various values of $\beta$ .

Figure 12.

Temperature profiles for the various values of $\beta$ .

Figure 13.

Velocity profiles for the various values of $\lambda$ .

Figure 14.

Temperature profiles for the various values of $\lambda$ .

Figure 15.

Velocity profiles for the various values of $S$ .

Figure 16.

Temperature profiles for the various values of $S$ .

Figure 17.

Temperature profiles for the various values of $M$ .

Figure 18.

Temperature profiles for the various values of $\mathrm{Rd}$ .

Figure 19.

Temperature profiles for the various values of $\mathrm{Bi}$ .

Figure 20.

Temperature profiles for the various values of $\mathrm{Pr}$ .

Conclusion

A numerical investigation is done to investigate the magnetohydrodynamic (MHD) mixed convection boundary layer flow and heat transfer characteristics of Fe₃O₄ and Ti₆Al₄V-Casson based nanofluid over an exponential stretching/shrinking surface. For this purpose, a pertinent flow model presented by Tiwari and Das ¹¹ is used to modify the governing equations. Similarity solutions of the governing equations are obtained through similarity transformation. The obtained equations are solved in the Maple software by shooting method. Triple solutions are found for different ranges of the used flow parameters. Based on numerical observations, following conclusions for the stable and physically feasible solution is made:

• The triple solutions are obtained at different ranges of suction parameter $\left(S\right)$ and the stretching/shrinking parameter $(\lambda )$ in relation with both assisting and opposing flow conditions.

• The rate of skin friction decreases for $\lambda >0$ and increases for $\lambda <0$ as the suction is increased for both Fe₃O₄ and Ti₆Al₄V-Casson based nanofluid flow.

• The heat transfer rate decreases by increasing in suction parameter for both Fe₃O₄ and Ti₆Al₄V-Casson nanofluid flow.

• Ti₆Al₄V-Casson nanofluid shows more inclined skin friction rate when compared to Fe₃O₄-Casson nanofluid while converse effect is observed in heat transfer rate.

• The velocity profile rises by increase in stretching/shrinking parameters and nanoparticles volumetric fractions.

• An increase in Casson and suction parameters decrease the corresponding velocity profiles. Also, their velocity boundary layer thicknesses reduce with increase in both parameters.

• The higher values of the Prandtl number, stretching/shrinking, and suction parameter which decrease the temperature profile.

• An increase in nanoparticles volumetric fraction, radiation and magnetic parameters as well as Biot number rise the temperature profiles and related thermal boundary layer thicknesses.