Effects of Stefan Blowing and Slip Conditions on Unsteady MHD Casson Nanofluid Flow Over an Unsteady Shrinking Sheet: Dual Solutions

Abstract

In this article, the magnetohydrodynamic (MHD) flow of Casson nanofluid with thermal radiation over an unsteady shrinking surface is investigated. The equation of momentum is derived from the Navier–Stokes model for non-Newtonian fluid where components of the viscous terms are symmetric. The effect of Stefan blowing with partial slip conditions of velocity, concentration, and temperature on the velocity, concentration, and temperature distributions is also taken into account. The modeled equations of partial differential equations (PDEs) are transformed into the equivalent boundary value problems (BVPs) of ordinary differential equations (ODEs) by employing similarity transformations. These similarity transformations can be obtained by using symmetry analysis. The resultant BVPs are reduced into initial value problems (IVPs) by using the shooting method and then solved by using the fourth-order Runge–Kutta (RK) technique. The numerical results reveal that dual solutions exist in some ranges of different physical parameters such as unsteadiness and suction/injection parameters. The thickness of the velocity boundary layer is enhanced in the second solution by increasing the magnetic and velocity slip factor effect in the boundary layer. Increment in the Prandtl number and Brownian motion parameter is caused by a reduction of the thickness of the thermal boundary layer and temperature. Moreover, stability analysis performed by employing the three-stage Lobatto IIIA formula in the BVP4C solver with the help of MATLAB software reveals that only the first solution is stable and physically realizable.

1. Introduction

In the past few decades, non-Newtonian fluids have attracted the interest of scientists, researchers, and mathematicians due to their significant applications in various industrial sectors. Many studies on non-Newtonian fluids have been carried out by considering several physical parameters. Tanveer et al. [1] found exact solutions for the non-Newtonian fluid known as Bingham nanofluid. Hayat et al. [2] examined non-Newtonian fluid under the influence of Brownian motion and the Bejan number. They concluded that the Bejan number and entropy rate display double behaviors in contradiction to the Eckert number. The Casson fluid model is one of the non-Newtonian fluid models in which yield stress characteristics are extensively studied. The fluids of this model fall within the dilatant liquid category. This model represents an inverse relationship between yield stress and shear stress. If the shear stress is less (more) than the yield stress, then the fluid behaves like a solid (fluid). Some common examples of this type are jam, sauce, nectar, soup, and organic product juices [3,4]. To the best of our knowledge, no such study on the unsteady flow of Casson nanofluid with the effect of Stefan blowing and partial slip conditions for multiple solutions has been reported in the literature.

Nowadays, nonlinear boundary value problems (BVPs) of fluid flows and their multiple solutions are important in engineering, physics, and mathematics due to their wide applications in engineering and scientific research. Therefore, it is very important not to miss any solution of nonlinear BVPs. Many numerical approaches fail to identify multiple solutions because the solutions may oscillate [5,6]. It is worth mentioning that the occurrence of multiple solutions of fluid flow over linear and nonlinear stretching surfaces is possible when the flow has a stagnation point or opposing flow [7]. On the contrary, flow over an exponential surface, which is quite different from the linear and nonlinear stretching surfaces, may have multiple solutions over a stretching surface even without cases of stagnation points and opposing flow [8]. In real situations, multiple solutions cannot be visualized in the boundary layer problems experimentally [9]. Therefore, mathematical analysis of the fluid flow model should be considered to detect the presence of multiple solutions. Terrill and Thomas [10] succeeded in discovering multiple solutions of laminar flow in a uniformly permeable pipe. In particular, their investigation uncovered that dual solutions exist outside the interval of estimations of suction Reynolds number 2.3 < R < 9.1 though no solution exists inside this range. Many other researchers [11,12,13] considered channel flow problems with uniformly permeable channels and reasoned that the existence and uniqueness theorem for their solutions is fulfilled for the scope of Reynolds number −∞ < R < ∞. Raithby [14] examined the problem of two-dimensional fluid flow in a permeable channel with suction and infusion numerically and found a second solution for R > 12. This investigation was restricted to various solutions for instances of suction R > 0. Lund et. al. [15] found multiple solutions of steady laminar Casson fluid with nanoparticles on a nonlinear stretching sheet with the assumption that an external magnetic field is applied on the flow. Their study revealed that multiple solutions exist only for the case of suction. Moreover, linear analysis of stability showed that only the first solution is physically reliable (stable). Meanwhile, Cu-C₆H₉N_aO₇ and Ag-C₆H₉N_aO₇ nanofluids were investigated by Lund et al. [16], who found that multiple solutions do not exist if the shrinking surface is impermeable. Recently, many authors [17,18,19,20,21,22] succeeded in obtaining multiple solutions of various fluid models under the influence of physical parameters. It was claimed that nonlinearity in the fluid model causes the existence of multiple solutions [23,24,25]. In light of the above discussion, it can be concluded that the existence of multiple solutions depends not only on nonlinearity in the equations but also on the values of different physical parameters.

Some applications of the magnetohydrodynamic (MHD) flow of Casson fluid can be seen in industrial sectors; therefore, the MHD characteristics of the flow need to be considered. Further, the Casson model shows distinct behavior when a magnetic field effect is imposed on it. The non-Newtonian problems of Casson fluid boundary layer flow over a shrinking surface have numerous applications in manufacturing processes and industry, particularly in the metal and plastic industries. There are various applications where a blowing effect is incorporated into the fluid flow problems. The effect of blowing on the MHD flow of a nanofluid on a shrinking surface has inspired numerous scholars to conduct further investigations. The Stefan blowing effect is different from mass blowing or injection because of transpiration [26]. Shahzad et al. [27] considered the MHD flow of Casson fluid with a thermal radiation effect. However, they found only one solution using the homotropy analysis method. Previously, multiple flow configurations involving the impact of blowing were explored by many researchers [22,28,29,30,31] for regular fluids. To the best of our knowledge, there have been no investigations focusing on the MHD flow of Casson nanofluids on a shrinking sheet with the Stefan blowing effect. Therefore, unsteady MHD flow of a Casson nanofluid with the combined effect of Stefan blowing and velocity slip conditions was taken into account in this study. Furthermore, thermal radiation and slip conditions were also considered.

2. Modeling and Simulation

The two-dimensional unsteady MHD flow of a Casson nanofluid on a shrinking surface, along with the effects of slip conditions and thermal radiation, was considered; these assumptions are presented in Figure 1.

The base (Casson) fluid with nanoparticles has fluid properties that are supposed to be not in thermal equilibrium. The isotropic equation of state is discussed thoroughly in Lund et al. [15,16] and Nakamura and Sawada [32], and can be expressed as

$${\tau }_{ij}=\{\begin{array}{l}\left({\mu }_B+\left(\frac{P_y}{\sqrt{2\pi }}\right)\right)2e_{ij},\pi >{\pi }_c\\ \left({\mu }_B+\left(\frac{P_y}{\sqrt{2{\pi }_c}}\right)\right)2e_{ij},\pi <{\pi }_c.\end{array}$$

(1)

It is also assumed that the flow is subjected to a transverse magnetic field of strength $B=\frac{B_0}{{\left(1-\epsilon t\right)}^{1/2}}$ where $B_0$ is the constant applied magnetic field. The effect of $B$ is applied perpendicular to the shrinking sheet (refer to Figure 1). Based on all assumptions, we have the following governing equations:

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

(2)

$$\frac{\partial \overline{u}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{v}}{\partial \overline{y}}=\vartheta \left(1+\frac{1}{\beta }\right)\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}-\frac{\sigma B^2\overline{u}}{\rho }$$

(3)

$$\frac{\partial \overline{T}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{T}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{T}}{\partial \overline{y}}=\alpha \frac{{\partial }^2\overline{T}}{\partial {\overline{y}}^2}+{\tau }_1\left[D_B\frac{\partial \overline{C}}{\partial \overline{y}}\frac{\partial \overline{T}}{\partial \overline{y}}+\frac{D_T}{T_{\infty }}{\left(\frac{\partial \overline{T}}{\partial \overline{y}}\right)}^2\right]-\frac{1}{{\left(\rho c_p\right)}_{nf}}{q}_r$$

(4)

$$\frac{\partial \overline{C}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{C}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{C}}{\partial \overline{y}}=D_B\frac{{\partial }^2\overline{C}}{\partial {\overline{y}}^2}+\frac{D_T}{T_{\infty }}\frac{{\partial }^2\overline{T}}{\partial {\overline{y}}^2}$$

(5)

subject to the boundary conditions

$$\begin{array}{l}\overline{v}=\frac{D_B}{1-{\overline{C}}_w}\frac{\partial \overline{C}}{\partial \overline{y}},\overline{u}={\overline{u}}_w+N_1\left(\overline{x},\overline{t}\right)\nu \left(1+\frac{1}{\beta }\right)\frac{\partial \overline{u}}{\partial \overline{y}},\overline{T}={\overline{T}}_w+D_1\left(\overline{x},\overline{t}\right)\frac{\partial \overline{T}}{\partial \overline{y}},\overline{C}={\overline{C}}_w+D_2\left(\overline{x},\overline{t}\right)\frac{\partial \overline{C}}{\partial \overline{y}}at\overline{y}=0\\ \overline{u}\to 0,\overline{T}\to {\overline{T}}_{\infty },\overline{C}\to {\overline{C}}_{\infty }\mathrm{as}\overline{y}\to \infty \end{array}$$

(6)

where $\beta ,\rho ,\mu ,k,\upsilon ,\mathrm{and}{\overline{u}}_w$ represent the Casson parameter, density, viscosity, thermal conductivity, and kinematic viscosity of the fluid, and the velocity of the plate, respectively. Further, $q_r=-\frac{4{\sigma }^{*}}{3K^{*}}\frac{\partial {\overline{T}}^4}{\partial \overline{y}},\alpha ,{\tau }_1\mathrm{and}{\left(\rho c_p\right)}_{nf}$ are, respectively, the radiation heat flux and thermal diffusivity of the Casson nanofluid, the ratio between the effective heat capacity of the nanoparticle material and the capacity of the fluid, and the heat capacitance in the nanofluid. $D_T$, $D_B,N_1\left(\overline{x},\overline{t}\right),D_1\left(\overline{x},\overline{t}\right),\mathrm{and}{D}_2\left(\overline{x},\overline{t}\right)$ are, respectively, the parameters of thermophoresis, mass diffusivity, velocity slip factor, thermal slip factor, and concentration slip factor. These parameters can be defined as

$$\{\begin{array}{c}\upsilon =\frac{\mu }{\rho },{\overline{u}}_w=-\frac{c\overline{x}}{1-\epsilon \overline{t}},N_1\left(\overline{x},\overline{t}\right)={\left(N_1\right)}_0\sqrt{1-\epsilon \overline{t}},\\ D_1\left(\overline{x},\overline{t}\right)={\left(D_1\right)}_0\sqrt{1-\epsilon \overline{t}},D_2\left(\overline{x},\overline{t}\right)={\left(D_2\right)}_0\sqrt{1-\epsilon \overline{t}}\\ \alpha =\frac{k_0}{{\left(\rho c_p\right)}_f},\beta =\frac{{\mu }_B\sqrt{2{\pi }_c}}{P_y},{\tau }_1=\frac{{\left(\rho c_p\right)}_p}{{\left(\rho c_p\right)}_f}.\end{array}$$

We use the following similarity transformations to get similarity solutions:

$$\eta ={\left(\frac{c}{\vartheta \left(1-\epsilon \overline{t}\right)}\right)}^{\frac{1}{2}}\overline{y},\overline{u}=\frac{c\overline{x}}{\left(1-\epsilon \overline{t}\right)}{f}^{\prime }\left(\eta \right),\overline{v}=-{\left(\frac{\vartheta c}{\left(1-\epsilon \overline{t}\right)}\right)}^{\frac{1}{2}}f\left(\eta \right),\theta \left(\eta \right)=\frac{\overline{T}-{\overline{T}}_{\infty }}{{\overline{T}}_w-{\overline{T}}_{\infty }},\varnothing \left(\eta \right)=\frac{\overline{C}-{\overline{C}}_{\infty }}{{\overline{C}}_w-{\overline{C}}_{\infty }}.$$

(7)

By substituting Equation (7) into Equations (3)–(6), the two-point BVPs below are obtained:

$$\left(1+\frac{1}{\beta }\right)f^{‴}-f^{\prime 2}+f{f}^{″}-A\left(0.5\eta f^{″}+f^{\prime }\right)-M{f}^{\prime }=0$$

(8)

$$\frac{\left(1+N_r\right)}{Pr}{\theta }^{″}+f{\theta }^{\prime }+N_t{{\theta }^{\prime }}^2+N_b{\theta }^{\prime }{\varnothing }^{\prime }-0.5A\eta {\theta }^{\prime }=0$$

(9)

$$\frac{1}{Sc}{\varnothing }^{″}+f{\varnothing }^{\prime }-0.5A\eta {\varnothing }^{\prime }+\frac{1}{Sc}\frac{N_t}{N_b}{\theta }^{″}=0$$

(10)

subject to newly transformed boundary conditions

$$\begin{gathered} f\left(0\right)=\frac{f_w}{Sc}{\varnothing }^{\prime }\left(0\right),f^{\prime }\left(0\right)=-1+\delta \left(1+\frac{1}{\beta }\right)f^{″}\left(0\right),\theta \left(0\right)=1+{\delta }_T{\theta }^{\prime }\left(0\right),\varnothing \left(0\right)=1+{\delta }_C{\varnothing }^{\prime }\left(0\right) \\ {f}^{\prime }\left(\eta \right)\to 0,\theta \left(\eta \right)\to 0,\varnothing \left(\eta \right)\to 0\eta \to \infty \end{gathered}$$

(11)

where $A$ is the unsteadiness parameter. It should be noted that $A<0$ indicates decelerating flow, while $A>0$ shows accelerating flow. $M,N_r,Pr,N_t,N_b,Sc,f_w,\delta ,{\delta }_T,\mathrm{and}{\delta }_C$ are the dimensionless magnetic parameter, thermophoresis parameter, thermal radiation parameter, Brownian motion parameter, suction/injection parameter, Schmidt number, thermal slip factor, velocity slip factor, and concentration slip factor, respectively. These dimensionless parameters are defined as

$$\{\begin{array}{c}M=\frac{\sigma B_0^2}{\rho c},A=\frac{\epsilon }{c},N_r=\frac{16{\sigma }^{*}{T}_{\infty }^3}{3K^{*}k},Pr=\frac{\nu }{\alpha },\\ N_t=\raisebox{1ex}{${\tau }_1\Delta T{D}_T$}\!\left/ \!\raisebox{-1ex}{$\vartheta {\overline{T}}_{\infty }$}\right.,N_b=\raisebox{1ex}{${\tau }_1{\overline{C}}_{\infty }{D}_B$}\!\left/ \!\raisebox{-1ex}{$\vartheta $}\right.,f_w=\raisebox{1ex}{${\overline{C}}_{\infty }$}\!\left/ \!\raisebox{-1ex}{$1-{\overline{C}}_w$}\right.\\ \delta ={\left(N_1\right)}_0\sqrt{c\vartheta },{\delta }_T={\left(D_1\right)}_0\sqrt{\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$\vartheta $}\right.},{\delta }_C={\left(D_2\right)}_0\sqrt{\raisebox{1ex}{$c$}\!\left/ \!\raisebox{-1ex}{$\vartheta $}\right.}.\end{array}$$

The physical quantities from an engineering point of view are the coefficient of skin friction ($C_{fx}$), local Nusselt ($N{u}_x$), and Sherwood number ($S{h}_x$), obtained from the following equations:

$$C_{fx}=\frac{{\tau }_w}{\rho {\overline{u}}_w{}^2},N{u}_x=\frac{\overline{x}{q}_1}{k\left(T_w-T_{\infty }\right)},S{h}_x=\frac{\overline{x}{q}_2}{k\left(C_w-C_{\infty }\right)}$$

(12)

where ${\tau }_w$, $q_1$, and $q_2$ are the shear stress, heat, and mass flux at the wall, respectively, represented as

$${\tau }_w=-\mu \left(1+\frac{1}{\beta }\right){\left(\frac{\partial \overline{u}}{\partial \overline{y}}\right)}_{y=0},q_1=-k{\left(\frac{\partial \overline{T}}{\partial \overline{y}}\right)}_{y=0}+{\left(q_r\right)}_{y=0},q_2=-k{\left(\frac{\partial \overline{C}}{\partial \overline{y}}\right)}_{y=0}.$$

(13)

Which leads to

$$C_{fx}\sqrt{R{e}_x}=\left(1+\frac{1}{\beta }\right)f^{″}\left(0\right),\frac{N{u}_x}{\sqrt{R{e}_x}}=-\left(1+N_r\right){\theta }^{\prime }\left(0\right),\frac{S{h}_x}{\sqrt{R{e}_x}}=-{\varnothing }^{\prime }\left(0\right)$$

(14)

where $R{e}_x=\raisebox{1ex}{${\overline{u}}_w\overline{x}$}\!\left/ \!\raisebox{-1ex}{$\nu $}\right.$ is the local Reynolds number.

3. Stability Analysis

When there exists more than one solution in any fluid flow problem, stability analysis must be performed in that study. It should be noted that we describe the first (second) solution instead of the upper (lower) branch throughout the whole manuscript. Lund et al. [33] and Dero et al. [34] stated in their papers that the only stable and physically possible solution is the first solution (upper branch), whereas the second solution (lower branch) is physically unrealizable and unstable. Lund et al. [35], Weidman and Sprague [36], Lund et al. [20], and others also considered this analysis in their studies. We need to introduce a new dimensionless time variable τ in order to perform a stability analysis of the solution where τ corresponds to the initial value problems (IVPs).

The new similarity transformation variables with τ and (7) can be written as

$$\begin{gathered} \overline{u}=\frac{c\overline{x}}{\left(1-\epsilon \overline{t}\right)}\frac{\partial f\left(\eta ,\tau \right)}{\partial \eta },\overline{v}=-{\left(\frac{\vartheta c}{\left(1-\epsilon \overline{t}\right)}\right)}^{\frac{1}{2}}f\left(\eta ,\tau \right),\theta \left(\eta ,\tau \right)=\frac{\overline{T}-{\overline{T}}_{\infty }}{{\overline{T}}_w-{\overline{T}}_{\infty }},\varnothing \left(\eta ,\tau \right)=\frac{\overline{C}-{\overline{C}}_{\infty }}{{\overline{C}}_w-{\overline{C}}_{\infty }} \\ \eta ={\left(\frac{c}{\vartheta \left(1-\epsilon \overline{t}\right)}\right)}^{\frac{1}{2}}\overline{y};\tau =\frac{c\overline{t}}{\left(1-\epsilon \overline{t}\right)}. \end{gathered}$$

(15)

By substituting Equation (15) into Equations (3)–(5), the following equations can be obtained:

$$\left(1+\frac{1}{\beta }\right)\frac{{\partial }^{3}f}{\partial {\eta }^3}-{\left(\frac{\partial f}{\partial \eta }\right)}^2+f\frac{{\partial }^{2}f}{\partial {\eta }^2}-A\left(0.5\eta \frac{{\partial }^{2}f}{\partial {\eta }^2}+\frac{\partial f}{\partial \eta }\right)-M\frac{\partial f}{\partial \eta }-\left(1+A\tau \right)\frac{{\partial }^{2}f}{\partial \tau \partial \eta }=0$$

(16)

$$\frac{\left(1+N_r\right)}{Pr}\frac{{\partial }^2\theta }{\partial {\eta }^2}+f\frac{\partial \theta }{\partial \eta }+N_t{\left(\frac{\partial \theta }{\partial \eta }\right)}^2+N_b\frac{\partial \varnothing }{\partial \eta }\frac{\partial \theta }{\partial \eta }-0.5A\eta \frac{\partial \theta }{\partial \eta }-\left(1+A\tau \right)\frac{\partial \theta }{\partial \tau }=0$$

(17)

$$\frac{1}{Sc}\frac{{\partial }^2\varnothing }{\partial {\eta }^2}+f\frac{\partial \varnothing }{\partial \eta }-0.5A\eta \frac{\partial \varnothing }{\partial \eta }+\frac{1}{Sc}\frac{N_t}{N_b}\frac{{\partial }^2\theta }{\partial {\eta }^2}-\left(1+A\tau \right)\frac{\partial \varnothing }{\partial \tau }=0$$

(18)

subject to the boundary conditions

$$\begin{gathered} f\left(0,\tau \right)=\frac{f_w}{Sc}\frac{\partial \varnothing \left(0,\tau \right)}{\partial \eta },\frac{\partial f\left(0,\tau \right)}{\partial \eta }=-1+\delta \frac{{\partial }^{2}f\left(0,\tau \right)}{\partial {\eta }^2},\theta \left(0,\tau \right)=1+{\delta }_T\frac{\partial \theta \left(0,\tau \right)}{\partial \eta },\varnothing =1+{\delta }_C\frac{\partial \varnothing \left(0,\tau \right)}{\partial \eta } \\ \frac{\partial f\eta \left(0,\tau \right)}{\partial \eta }=\theta \left(\eta ,\tau \right)=\varnothing \left(\eta ,\tau \right)=0\hspace{1em}as\hspace{1em}\eta \to \infty . \end{gathered}$$

(19)

To check the stability of steady flow solutions where $f\left(\eta \right)=f_0\left(\eta \right),\theta \left(\eta \right)={\theta }_0\left(\eta \right)$, and $\varnothing \left(\eta \right)={\varnothing }_0\left(\eta \right)$ satisfying the boundary value problem (8)–(11), we write

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

(20)

where $F\left(\eta \right)$, $G\left(\eta \right)$, and $H\left(\eta \right)$ are small relative to $f_0\left(\eta \right)$, ${\theta }_0\left(\eta \right)$, and ${\varnothing }_0\left(\eta \right)$, respectively. Further, γ is the unknown eigenvalues. An infinite set of eigenvalues is obtained by solving the eigenvalue problem (21)–(23). From that set, we have to choose the smallest values of γ. Using Equation (20) in Equations (16)–(19) leads to

$$\left(1+\frac{1}{\beta }\right)F_0^{‴}-2f_0^{\prime }{F}_0^{\prime }+f_0{F}_0^{″}+F_0{f}_0^{″}-A\left(0.5\eta F_0^{″}+F_0^{\prime }\right)-M{F}_0^{\prime }+\mathsf{\gamma }{F}_0^{\prime }=0$$

(21)

$$\frac{\left(1+N_r\right)}{Pr}{G}_0^{″}+f_0{G}_0^{\prime }+F_0{\theta }_0^{\prime }+2N_t{\theta }_0^{\prime }{G}_0^{\prime }+N_b\left({\varnothing }_0^{\prime }{G}_0^{\prime }+H_0^{\prime }{\theta }_0^{\prime }\right)-0.5\eta A{G}_0^{\prime }+\mathsf{\gamma }{G}_0=0$$

(22)

$$H_0^{″}+Sc\left(f_0{\varnothing }_0^{\prime }+F_0{H}_0^{\prime }\right)+\frac{N_t}{N_b}{G}_0^{″}-0.5\eta ASc{H}_0^{\prime }+Sc\mathsf{\gamma }{H}_0=0$$

(23)

subject to boundary conditions

$$\begin{gathered} F_0\left(0\right)=\frac{f_w}{Sc}{H}_0\left(\eta \right),F_0^{\prime }\left(0\right)=\delta \left(1+\frac{1}{\beta }\right)F_0^{″}\left(0\right),G_0\left(0\right)={\delta }_T{G}_0^{\prime }\left(0\right),H_0\left(0\right)={\delta }_C{H}_0^{\prime }\left(0\right), \\ {F}_0^{\prime }\left(\eta \right)\to 0,G_0\left(\eta \right)\to 0,H_0\left(\eta \right)\to 0,as\eta \to \infty . \end{gathered}$$

(24)

According to Lund et al. [37], in order to assess the stability of Equations (21)–(24) we need to relax one boundary condition on $F_0^{\prime }\left(\eta \right),G_0\left(\eta \right)$, or $H_0\left(\eta \right)$. It should be noted that we relaxed $F_0^{\prime }\left(\eta \right)\to 0as\eta \to \infty$ to $F_0^{″}\left(0\right)=1$ in this problem.

4. Results and Discussion

In this study, the unsteady MHD flow of a Casson nanofluid on a shrinking sheet under the influence of thermal radiation and slip conditions was examined. The graphical results of the velocity, concentration, and temperature distributions for magnetic parameter $M$, thermal radiation parameter $N_r$, thermophoresis parameter $N_t,$ Brownian motion $N_b$, Schmidt number $Sc$, parameter of suction/injection $f_w$, velocity slip factor $\delta$, thermal slip factor ${\delta }_T$, and concentration slip factor ${\delta }_C$ were taken into account.

4.1. Analysis of Skin Friction, Temperature, and Concentration Rates

A comparison of the numerical values of the coefficient of skin friction obtained by fixing the values of $M=0.5,N_r=0.5,Pr=1,N_t=0.5,N_b=0.3,Sc=1,\delta =0.1,{\delta }_T=0.1,\mathrm{and}{\delta }_C=0.1$ is presented in

Table 1. Comparison of the shooting method with the three-stage Lobatto IIIA formula for values of $f^{″}\left(0\right)$ when $M=0.5,N_r=0.5,Pr=1,N_t=0.5,N_b=0.3,Sc=1,\delta =0.1,{\delta }_T=0.1,\mathrm{and}{\delta }_C=0.1$.

${\mathit{f}}_{\mathit{w}}$	$\mathit{\beta }$	$\mathit{A}$	Shooting Method		Three-Stage Lobatto	IIIA Formula
			1st Solution	2nd Solution	1st Solution	2nd Solution
−5	$\infty$	−1	5.2069632	0.6228652	5.2069632	0.6228652
−5.1			5.8082723	0.5521204	5.8082723	0.5521204
−5.2	$\infty$		6.1370006	0.4934203	6.1370006	0.4934203
	1.5		3.6769605	1.7569470	3.6769605	1.7569470
	2.5	−1	4.3796524	1.8324565	4.3796524	1.8324565
		−0.5	4.3574468	2.6719383	4.3574468	2.6719383
		0	4.3317923	2.9915082	4.3317923	2.9915082
		0.5	4.3013107	3.2212830	4.3013107	3.2212830
		1	4.2634822	3.4128383	4.2634822	3.4128383

. From this table, we found an excellent agreement with the numerical values of the skin friction coefficient computed using two different numerical approaches. The collective effects of suction/injection parameter $f_w$ and Casson parameter $\beta$ on the skin friction $f^{″}\left(0\right)$ are shown in Figure 2. It was observed from this figure that there exist dual solutions for Casson parameters $\beta =1.5$ and $\beta =2.5$ at $f_w=-4.9909$ and $f_w=-4.98645,$ respectively. Without any doubt, we can say that $f_w=-4.9909$ and $f_w=-4.98645$ are our critical points where multiple solutions exist. In the same vein, we concluded from this profile that the skin friction coefficient magnitude increases monotonically for both solutions (first solution and second solution).

Figure 3 exhibits the skin friction coefficient with variation in unsteady parameter $A$ for suction parameters $f_w=-5.2,-5.3,-5.4$. The skin friction increases gradually for the first solution when the values of the unsteady parameter and suction parameter are increased. However, a completely contradictory trend was observed from the same profile for the second solution. It can be concluded from Figure 3, Figure 4 and Figure 5 that there exist two ranges of solutions, namely, no solution and dual solutions. Dual solutions exist when $A\le A_i$ where $i=1,2,3$. When $A>A_i$, a solution does not exist over the shrinking surface. Figure 4 demonstrates the heat transfer rates for various values of unsteady parameter and suction parameter. This figure shows that the heat transfer rate increases monotonically with increasing suction and unsteady parameters. This is because unsteadiness leads the fluid particles to migrate from one place to another with a collisional effect which produces heat. Therefore, heat transfer increases (decreases) gradually as the value of the unsteady parameter increases in the first (second) solution. Figure 5 shows the effect of the unsteady parameter on the concentration rate. It can be observed from this profile that the concentration transfer rate increased (decreased) with increasing unsteady parameter and suction parameter for the first (second) solution.

4.2. Analysis of Velocity Profiles

The velocity profiles under the variation of various physical parameters are plotted in Figure 6, Figure 7 and Figure 8. Figure 6 demonstrates the velocity profile with variation of the Casson magnetic parameters. It can be observed that the velocity profile for both solutions decreases with enhanced magnetic field strength for both cases, Newtonian ($\beta \to \infty$) and non-Newtonian ($\beta =2.5$). In practical terms, it can be said that an increase in the magnetic parameter generates more Lorentz force, which then reduces the velocity of the fluid due to drag force strength. Consequently, the velocity profile and thickness of the momentum boundary layer decline as the magnetic parameter increases. The impact of the velocity slip parameter on the velocity profile is depicted in Figure 7. From this figure, we came to understand that the velocity profile inclines (declines) for the first (second) solution as the velocity slip parameter increases. Figure 8 illustrates the effect of the unsteady parameter on the profile of velocity for Newtonian ($\beta \to \infty$) and non-Newtonian ($\beta =2.5$) cases. This profile shows that the momentum boundary layer thickness increases gradually with increasing value of the unsteady parameter for both Newtonian and non-Newtonian cases. Therefore, the velocity profile increases for both solutions.

4.3. Analysis of Temperature Profiles

Figure 9 presents the behavior of the temperature profile for various values of Prandtl number in the presence of a thermal radiation effect. It was observed from this profile that the radiation parameter tends to increase the temperature of the nanofluid as the numerical value of the radiation parameter is enhanced. Physically, we can say that the strength of the radiation parameter tends to develop the thermal boundary layer. However, the opposite trend was noticed for higher values of the Prandtl number. According to Khan et al. [38], “an increase in Pr accompanies with weaker thermal diffusivity and restricts the heat from flowing deeper into the nanofluid, so thermal boundary layer becomes decreased with an augmentation of Pr”.

The impact of the thermal slip parameter on the temperature profile is depicted in Figure 10. It can be observed that the thickness of the thermal boundary layer decreases (increases) for the first (second) solution. As a result, the temperature profile decreases for the first solution, and the opposite trend was observed for the second solution. The temperature profile for both solutions is gradually reduced as the Brownian motion parameter $N_b$ increases (see Figure 11). The effect of thermophoresis $N_t$ on the temperature profile is shown in Figure 12. A force, known as thermophoresis force, pushes the nanoparticles towards the ambient flow, and as a result, the thermal boundary layer becomes thicker. It is worth mentioning that the temperature gradient becomes smaller due to the huge number of nanoparticles in the fluid flow and thus increases the thermal boundary layer thickness and nanoparticle volume fraction.

4.4. Analysis of Concentration Profiles and Stability

Figure 13 and Figure 14 show the effect of $N_b$ and $N_t$ on the concentration profiles. In these figures, it can be easily seen that thermophoresis and Brownian motion parameters have an inverse relation to each other with regard to the distributions of concentration. The thickness of the concentration layer becomes thinner as the Brownian motion impact is increased. Physically, this is due to the fact that Brownian motion creates more random movement of nanoparticles inside the fluid flow, and as a result, the concentration of nanoparticles decreases. However, the opposite trend was noticed in Figure 14.

The values of the smallest eigenvalues $\gamma$ for various values of $f_w$ are presented in Figure 15. Linear stability analysis was adopted to plot this graph, which is divided into two regions. The lower half of the region is the area where the negative eigenvalues are plotted, which refers to an initial growth of disturbance, and flow is therefore unstable in this mode. However, in the upper half, the positive eigenvalues show an initial decay of the disturbance, and therefore, the flow is stable. In the same vein, the smallest eigenvalues tend to zero for the upper and lower branches of the solutions as $f_w\to -4.98645$. Physically, we can say that $f_w=-4.98645$ is the critical point where the system has multiple solutions. Moreover, from this graphical representation, we can conclude that the first solution is physically reliable and stable. On the other hand, the second solution is physically unreliable and unstable.

5. Conclusions

In this study, we considered the unsteady flow of electrically conductive Casson nanofluid in the presence of thermal radiation over an unsteady shrinking surface. The effect of Stefan blowing with partial slip conditions of velocity, temperature, and concentration on the velocity, temperature, and concentration distributions was examined. The modeled PDEs were transformed into ODEs in the form of BVPs by employing similarity transformations. The resultant BVPs were reduced to IVPs using the shooting method and then solved using the fourth-order Runge–Kutta (RK) technique. The main pointwise conclusions that can be drawn from this investigation are as follows:

• Dual solutions exist when $A\le A_i$ where $i=1,2,3$ and $f_w\le f_{wc}$ where $c=1,2$.
• The temperature and thickness of the thermal boundary layer are reduced when the Prandtl number and Brownian motion parameter are increased.
• The velocity boundary layer becomes thicker in the second solution when the magnetic and velocity slip factor effect is increased.
• Beyond the critical point, there is a range with no solution.
• The results of stability analysis revealed that there exists initial decay (growth) of disturbances for the first (second) solution.
• The more physical realizable solution is the first solution.