Triple Solutions and Stability Analysis of Micropolar Fluid Flow on an Exponentially Shrinking Surface

Abstract

In this article, we reconsidered the problem of Aurangzaib et al., and reproduced the results for triple solutions. The system of governing equations has been transformed into the system of non-linear ordinary differential equations (ODEs) by using exponential similarity transformation. The system of ODEs is reduced to initial value problems (IVPs) by employing the shooting method before solving IVPs by the Runge Kutta method. The results reveal that there are ranges of multiple solutions, triple solutions, and a single solution. However, Aurangzaib et al., only found dual solutions. The effect of the micropolar parameter, suction parameter, and Prandtl number on velocity, angular velocity, and temperature profiles have been taken into account. Stability analysis of triple solutions is performed and found that a physically possible stable solution is the first one, while all leftover solutions are not stable and cannot be experimentally seen.

1. Introduction

Generally, the investigation of the non-Newtonian fluid flow in two-dimensional problems is a hard task when multiple solutions are attempted to find because equations are related to high nonlinear terms. In spite of these challenges, the researchers are making efforts to tackle these problems for multiple solutions due to their wide range of applications in various science and industrial fields. The category of non-Newtonian fluids which pacts with the suspended micro-rotational particles, is known as the micropolar fluid. Eringen introduced the theory of micropolar [1,2]. He explained the impacts of couple stresses and local rotational inertia that cannot be described by the standard equations of Navier–Stokes. The micropolar equations are mathematically described for the theory of porous media and the theory of lubrication in the books by Lukaszewicz [3] and Eringen [4]. There are various applications of the micropolar fluids, for example, liquid crystals, particle suspension, animal blood, lubrication, turbulent shear flows, and paints. Lok et al., [5] considered the stagnation point flow of micropolar nanofluid and then succeeded to find the dual solution. It is also stated that only the stable solution is the upper branch by performing an analysis of stability. Sheremet et al., [6] examined micropolar fluid with a convectional effect in a wavy triangular cavity. Bhattacharyya et al., [7] studied the micropolar fluid with thermal radiation effect and discovered that there exist dual solutions, while they did not perform the stability analysis. The dual solution has been discovered in the opposing flow case for the micropolar fluid [8], whereas stability analysis has not been performed in their examination of the multiple solutions. Ramzan et al., [9] investigated the nanofluid of micropolar non-Newtonian kind of fluid on the stretching sheet and claimed that velocity and angular velocity have an inverse relation with slip parameter. Turkyilmazoglu [10] analytically examined the Magnetohydrodynamic (MHD) flow of micropolar fluid. Shah et al., [11] examined micropolar nanofluid with the effect of the Casson parameter in the channel and stated that the thermal boundary layer becomes thicker for the higher values of the Brownian motion parameter. Some papers have recently shown some development on the micropolar fluid [12,13,14,15,16].

Nowadays, multiple solutions of fluid flow have received extensive consideration. Generally, many researchers stated in their studies that non-uniqueness of solutions, depending on the non-linearity in the fluid flow equations. Jawad et al. [17] stated that, due to the nonlinearity in fluid flow problems, there is more than one solution; however, seeking all possible solutions challenges researchers. Further, several of the approaches are ineffective because solutions are close together. Rohni [18] reported in her Ph.D. thesis that “the multiplicity of solutions in fluid dynamics and heat transfer is important to be able to be computed, since solutions arising from bifurcations often interact with one another producing otherwise inexplicable phenomena”. Moreover, she claimed that, in order to notice all possible solutions of any fluid, the flow problem is still a challenge. Mishra and DebRoy [19] described that multiple solutions have many important applications when these are related to heat transfer, because the final qualities and structure of many products of material processing in the industries can be improved by the concept of multiple solutions. Lund et al., [20] considered nanofluid, with the base fluid as sodium alginate and solid particles of copper and silver, over the shrinking sheet and obtained dual solutions and reported that the existence of multiple solutions depends upon how the researchers set the ranges of different applied physical parameters. Furthermore, they specified that the dual solution depended on the range of the suction parameter. Khashi’ie et al. [21] highlighted the importance of the multiple solutions and the stability analysis during the examination of the mixed convection flow of the micropolar fluid. They stated “if the problem has non-unique solutions but the researchers manage to find one solution only, there is a probability that the solution is the lower branch solution (unstable/not real). This will lead to the misinterpretation of the flow and heat transfer characteristics”. These papers include some significant studies on multiple solutions [22,23,24]. Here, it is worth highlighting that, when multiple solutions appeared in any fluid flow problem, the study of stability analysis should be considered. The first step in performing the stability analysis is, according to Rana et al. [25], to change the governing equation to unsteady form by adding a new time dependent variable for similarity. Lund et al., [26] found the dual solution of the Williamson fluid and claimed that there existed an infinite number of the values of the smallest eigenvalue by performing the stability of solutions. Waini et al., [27] reported that the stability of any solution dependent upon the sign of the value of the smallest eigenvalue. If the value of eigenvalue is positive, it indicated the stable solution, while the negative value of the smallest eigenvalue indicated the unstable solution [28].

The pointwise prime objectives of this study are, as follows:

• Reconsidering the problem of Aurangzaib et al. [29].
• To find all possible multiple solutions.
• To perform the stability analysis that has not been considered by the Aurangzaib et al. [29].
• Indicating the stable and unstable solutions by doing stability analysis, which cannot be experientially seen, due to that mathematical analysis is necessary.

The authors tried their best to find all possible solutions due to the importance of multiple solutions and lack of available literature in which triple solutions of any fluid flow problem were noticed. This paper is divided into six sections; Section 1 is for the brief introduction of micropolar fluid and multiple solutions. Mathematical formulation, derivation of stability analysis, and methodology are kept in Section 2, Section 3, and Section 4 respectively. Section 5 is for the result and discussion and Section 6 is for the conclusion.

2. Mathematical Formulation

An incompressible laminar boundary layer two-dimensional flow of the micropolar fluid over the exponentially shrinking sheet has been considered. The corresponding velocities of x and y-axes are u and v. The shrinking velocity is assumed to be $u_w\left(x\right)=-U_w{e}^{\frac{2x}{\ell }}$. The temperature of the sheet is taken to be $T_{ w}\left(x\right)=T_{\infty }+T_0{e}^{\frac{x}{2\ell }}$, as shown in Figure 1. The N = N(x,y) is supposed as the angular velocity. The respective boundary layer movement equation, along with micro rotations and the heat transfer equations can be expressed as vectors in accordance with the abovementioned assumptions.

$$\nabla ·\mathit{V}=0$$

(1)

$$\rho \frac{d\mathit{V}}{dt}=-\nabla \mathit{P}+\left(\mu +\kappa \right){\nabla }^2\mathit{V}+\kappa \left(\nabla \times \mathit{N}\right)$$

(2)

$$\rho j\frac{d\mathit{N}}{dt}=\gamma {\nabla }^2\mathit{N}-\kappa \left(2\mathit{N}-\nabla \times \mathit{V}\right)$$

(3)

$$\rho c_p\frac{dT}{dt}=k{\nabla }^{2}T$$

(4)

in which velocity vector is $\mathit{V}\equiv \left[u\left(x,y\right),v\left(x,y\right),0\right]$, the micro-rotation vector is $\mathit{N}$, ρ stands for fluid density, µ for viscosity coefficient, $\kappa$ for vertex viscosity, j is the density of micro-rotation, and γ stands for micropolar constant. We get following boundary layer equations according to the scale analysis.

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

(5)

$$u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}=\left(\vartheta +\frac{\kappa }{\rho }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\kappa }{\rho }\frac{\partial N}{\partial y}$$

(6)

$$u\frac{\partial N}{\partial x}+v\frac{\partial N}{\partial y}=\frac{\gamma }{\rho j}\frac{{\partial }^{2}N}{\partial y^2}-\frac{\kappa }{\rho j}\left(2N+\frac{\partial u}{\partial y}\right)$$

(7)

$$u \frac{\partial T}{\partial x}+v \frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2}$$

(8)

Subject to these boundary conditions

$$\begin{gathered} v=v_w\left(x\right); u=u_w\left(x\right); N=-m\frac{\partial u}{\partial y}; T=T_{ w}\left(x\right); at y=0 \\ u\to 0; N\to 0; T\to T_{\infty } as y\to \infty \end{gathered}$$

(9)

Now, we look for similarity transformation variables in order to transform Equations (6)–(8) with boundary conditions (9)

$$\begin{gathered} u=U_w{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$l$}\right.}{f}^{\prime }\left(\eta \right);v=-\sqrt{\frac{\vartheta U_w}{2l}}{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.}\left(f\left(\eta \right)+\eta f^{\prime }\left(\eta \right)\right);N=U_w{e}^{\raisebox{1ex}{$3x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.}\sqrt{\frac{U_w}{2\vartheta l}}h\left(\eta \right); \\ \theta \left(\eta \right)=\raisebox{1ex}{$\left(T-T_{\infty }\right)$}\!\left/ \!\raisebox{-1ex}{$\left(T_w-T_{\infty }\right)$}\right.; \eta =y\sqrt{\frac{U_w}{2\vartheta l}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2l$}\right.} \end{gathered}$$

(10)

By applying Equation (10) in Equations (6)–(9), we have the following system of similarity transformed ordinary differential equations

$$\left(1+K\right)f^{‴}+f{f}^{″}-2f^{’2}+K{h}^{\prime }=0$$

(11)

$$\left(1+\frac{K}{2}\right)h^{″}+f{h}^{\prime }-3f^{\prime }h-K\left(2h+f^{″}\right)=0$$

(12)

$$\frac{1}{Pr}{\theta }^{″}+f{\theta }^{\prime }-f^{\prime }\theta =0$$

(13)

Subject to boundary conditions

$$\begin{array}{c}f\left(0\right)=S; f^{\prime }\left(0\right)=-1; h\left(0\right)=-n{f}^{″}\left(0\right); \theta \left(0\right)=1\\ f^{\prime }\left(\eta \right)\to 0; h\left(\eta \right)\to 0; \theta \left(\eta \right)\to 0; as \eta \to \infty \end{array}$$

(14)

where prime denotes the differentiation with respect to $\eta$, the micropolar material parameter is $K=\frac{\kappa }{\mu }$, Prandtl number is $Pr=\frac{\vartheta }{\alpha }$, and suction is $S=-\frac{v_w}{\sqrt{\vartheta U_w/2l}}$.

The physical quantities of interest include skin friction, the stress of local couples, and the local number of Nusselt, which are described as

$$C_f=\frac{{\left[\left(\mu +\kappa \right)\frac{\partial u}{\partial y}+\kappa N\right]}_{y=0}}{\rho u_w^2};M_x=\frac{-\gamma {\left[\frac{\partial N}{\partial y}\right]}_{y=0}}{\rho x{u}_w^2}; N_u=\frac{-x{\left[\frac{\partial T}{\partial y}\right]}_{y=0}}{\left(T_w-T_{\infty }\right)}$$

(15)

By applying similarity transformation variables (10) in Equation (15), we have

$$\begin{gathered} C_f{\left(R{e}_x\right)}^{\frac{1}{2}}\sqrt{2l/x}=\left(1+\left(1-m\right)K\right)f^{″}\left(0\right),M_{x}R{e}_x=\left(1+\frac{K}{2}\right)h^{\prime }\left(0\right), \\ {N}_u{\left(R{e}_x\right)}^{-\frac{1}{2}}\sqrt{2l/x}=-{\theta }^{\prime }\left(0\right) \end{gathered}$$

(16)

where $R{e}_x=x{u}_w/\vartheta$ is the local Reynolds number.

3. Stability Analysis

According to Nasir et al., [28] and Rana et al., [25], we need to introduce the unsteady form of Equations (6)–(8) in order to perform stability test,

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}=\left(\vartheta +\frac{\kappa }{\rho }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\kappa }{\rho }\frac{\partial N}{\partial y}$$

(17)

$$\frac{\partial N}{\partial t}+u\frac{\partial N}{\partial x}+v\frac{\partial N}{\partial y}=\frac{\gamma }{\rho j}\frac{{\partial }^{2}N}{\partial y^2}-\frac{\kappa }{\rho j}\left(2N+\frac{\partial u}{\partial y}\right)$$

(18)

$$\frac{\partial T}{\partial t}+u \frac{\partial T}{\partial x}+v \frac{\partial T}{\partial y}=\frac{k^{*}}{\rho c_p}\frac{{\partial }^{2}T}{\partial y^2}$$

(19)

where time is denoted by t. The similarity transformations are written according to Roşca and Pop [30], as,

$$\begin{gathered} \psi =\sqrt{2\vartheta La}{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2L$}\right.}f\left(\eta ,\tau \right); N=a{e}^{\raisebox{1ex}{$3x$}\!\left/ \!\raisebox{-1ex}{$2L$}\right.}\sqrt{\frac{a}{2\vartheta L}}h\left(\eta ,\tau \right); \\ \theta \left(\eta ,\tau \right)=\raisebox{1ex}{$\left(T-T_{\infty }\right)$}\!\left/ \!\raisebox{-1ex}{$\left(T_w-T_{\infty }\right)$}\right.;\eta =y\sqrt{\frac{a}{2\vartheta L}} e^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$2L$}\right.} ; \tau =\frac{a}{2L}{e}^{\raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}.t \end{gathered}$$

(20)

By applying these similarity transformations, we reduced our Equations (17)–(19) into the following form

$$\left(1+K\right)\frac{{\partial }^{3}f}{\partial {\eta }^3}+f\frac{{\partial }^{2}f}{\partial {\eta }^2}-2{\left(\frac{\partial f}{\partial \eta }\right)}^2+K\frac{\partial h}{\partial \eta }-\frac{{\partial }^{2}f}{\partial \tau \partial \eta }=2\tau \left(\frac{\partial f}{\partial \eta }\frac{{\partial }^{2}f}{\partial \eta \partial x}-\frac{{\partial }^{2}f}{\partial {\eta }^2}\frac{\partial f}{\partial x}\right)$$

(21)

$$\left(1+\frac{K}{2}\right) \frac{{\partial }^{2}h}{\partial {\eta }^2}+f\frac{\partial h}{\partial \eta }-3h\frac{\partial f}{\partial \eta }-2Kh-K\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{\partial h}{\partial \tau }=2\tau \left(\frac{\partial f}{\partial \eta }\frac{\partial h}{\partial x}-\frac{\partial h}{\partial \eta }\frac{\partial f}{\partial x}\right)$$

(22)

$$\frac{1}{Pr}\frac{{\partial }^2\theta }{\partial {\eta }^2}+f\frac{\partial \theta }{\partial \eta }-\theta \frac{\partial f}{\partial \eta }-\frac{\partial \theta }{\partial \tau }=2\tau \left(\frac{\partial f}{\partial \eta }\frac{\partial \theta }{\partial x}-\frac{\partial \theta }{\partial \eta }\frac{\partial f}{\partial x}\right)$$

(23)

with the following boundary conditions

$$\begin{gathered} f\left(0,\tau \right)=S; \frac{\partial f\left(0, \tau \right)}{\partial \eta }=-1; h\left(0,\tau \right)=-m\frac{{\partial }^{2}f\left(0, \tau \right)}{\partial {\eta }^2}; \theta \left(0,\tau \right)=1 \\ \frac{\partial f\left(\eta , \tau \right)}{\partial \eta }\to 0; h\left(\eta ,\tau \right)\to 0; \theta \left(\eta ,\tau \right)\to 0 as \eta \to \infty \end{gathered}$$

(24)

and then we applied perturbed with the disturbance (see Roşca and Pop [30]) in Equations (21)–(24) with the following functions

$$\begin{gathered} f\left(\eta ,\tau \right)=f_0\left(\eta \right)+e^{-\epsilon \tau }F\left(\eta ,\tau \right) \\ h\left(\eta ,\tau \right)=h_0\left(\eta \right)+e^{-\epsilon \tau }G\left(\eta ,\tau \right) \\ \theta \left(\eta ,\tau \right)={\theta }_0\left(\eta \right)+e^{-\epsilon \tau }H\left(\eta ,\tau \right) \end{gathered}$$

(25)

where $F\left(\eta ,\tau \right)$, $G\left(\eta ,\tau \right)$, and $H\left(\eta ,\tau \right)$are small relative to $f_0\left(\eta \right)$, $h_0\left(\eta \right)$, and ${\theta }_0\left(\eta \right)$, respectively. Further, $\epsilon$ is the unknown eigenvalue. By substituting Equation (25) in Equations (21)–(23) by keeping $\tau =0$, we have following the linearized eigenvalue problems

$$\left(1+K\right)F_0^{’’’}+f_0{F}_0^{’’}+F_0{f}_0^{’’}-4f_0^{’}{F}_0^{’}+K{G}_0^{’}+\epsilon F_0^{’}=0$$

(26)

$$\left(1+\frac{K}{2}\right)G_0^{’’}+f_0{G}_0^{’}+F_0{h}_0^{’}-3h_0{F}_0^{’}-3h_0{F}_0^{’}-2K{G}_0-K{F}_0^{’’}+\epsilon G_0=0$$

(27)

$$\frac{1}{Pr}{H}_0^{’’}+f_0{H}_0^{’}+F_0{\theta }_0^{’}-f_0^{’}{H}_0-F_0^{’}{\theta }_0+\epsilon H_0=0$$

(28)

With the boundary conditions

$$\begin{array}{c}{F}_0\left(0\right)=0, F_0^{’}\left(0\right)=0, G_0\left(0\right)=-n{F}_0^{’’}\left(0\right), H_0\left(0\right)=0\\ F_0^{’}\left(\eta \right)\to 0, G_0\left(\eta \right)\to 0, H_0\left(\eta \right)\to 0, as \eta \to \infty \end{array}$$

(29)

We have to solve above linearized Equations (10)–(13) with new relax boundary conditions in order to find the values of smallest eigenvalue. In this particular problem, we have relaxed $H_0\left(\eta \right)\to 0$ as $\eta \to \infty$ into $H_0^{\prime }\left(0\right)=1$, see [31,32,33].

4. Numerical Methods

The governing ODEs are highly non-linear and, therefore, we adopt the numerical approach in order to solve Equations (11)–(14) and Equations (26)–(19). In this study, two methods have been employed, namely shooting method and Three-stage Lobatto III-A formula, which were used in many research articles of same authors previously [refer to 17,20,23,26]. The descriptions regarding these methods are explained below.

4.1. Shooting Method

The shooting technique along with the Runge Kutta method of the fourth order is employed in order to obtain the numerical solutions of Equations (11)–(13) subject to the boundary conditions. Shooting method helps to reduce the third order ODEs (11)–(13) into the first-order ODEs, such that

$$p=f’,q=p’, r=h’; q’=\frac{1}{1+K}\left\{2p^2-fq-Kr\right\}$$

(30)

$$r=h; r’=\frac{2}{2+K}\left\{3ph-fr+K\left(2h+q\right)\right\}$$

(31)

$$s=\theta ’; s’=Pr\left\{p\theta -fs\right\}$$

(32)

with conditions

$$f\left(0\right)=S; p\left(0\right)=-1; q\left(0\right)={\alpha }_1; h\left(0\right)=-n{\alpha }_1;r\left(0\right)={\alpha }_2; \theta \left(0\right)=1;s\left(0\right)={\alpha }_3$$

where ${\alpha }_1$, ${\alpha }_2,$ and ${\alpha }_3$ are called as unknown initial conditions. These three missing values ${\alpha }_1$, ${\alpha }_2,$ and ${\alpha }_3$ have to be obtained by using different shoots; this process of shoots will be continue until the profiles of the $f^{\prime }\left(\eta \right)\to 0; h\left(\eta \right)\to 0;\mathrm{and} \theta \left(\eta \right)\to 0$are satisfied the boundary condition$\eta \to \infty$. Maple (18) software has been used to convert the system of the third order ODEs into the system of the first order ODEs, for this process shootlib function is built-in Maple. Using RK method solves the system of the first order ODEs. Further, a detailed discussion about the shooting method with Maple software can be seen in the paper of Meade et al. [34].

4.2. Three-Stage Lobatto III-A Formula

Three-stage Lobatto III-A formula is built in BVP4C function with aid of $C^1$ piece-wise cubic polynomial in the finite difference code. According to Lund et al., [35] and Raza et al., [36], “this collocation polynomial and formula offers a $C^1$ continuous solution in which mesh error control and selection are created on the residual of the continuous solution. The tolerance of relative error is fixed ${10}^{-5}$ for the current problem. The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns solution, called as sol.y., as a construction. In any case, values of the solution are gotten from the array named sol.y relating to the field sol.x”. In addition, Figure 2 explains the algorithm of the method for stability analysis of the solutions.

5. Result and Discussion

Aurangzaib et al., [29] solved these Equations (11)–(14) and found dual solutions, which was the main contribution of authors. Aurangzaib et al., [29] have given some strong statements during the investigation of critical points. Before going to detail, we check the accuracy of our method (shooting method) by comparing our results with previously published literature in

Table 1. The comparison of values of heat transfer rate for different values of Pr.

Pr	M	Ishak [33]	Pramanik [37]	Raju et al. [38]	Present Results
1	0	0.9548	0.9547	0.954734	0.954955
2	0	1.4715	1.4714	1.471426	1.471421
3	0	1.8691	1.8691	1.869134	1.869044
10	0	3.6603	3.6603	3.660312	3.660354

and found excellent agreement with them. Moreover, we also qualitatively compared our results with Aurangzaib et al., [29] in Figure 3 and found in the good agreement. This gives trust in our numerical calculation and urges us to additionally contemplate this problem. Aurangzaib et al., [29] found dual solutions and stated in the result and discussion section that “the present study shows that for K = 0.1, i.e., for micropolar fluid, the similarity solutions exist when $\mathrm{S}\ge 2.3231$ and no similarity solution exists for $\mathrm{S}<2.3231$”. Firstly, we would like to clarify that there exist triple solutions not dual solutions; secondly, there is a range of multiple solutions and single solution. From Figure 4, a conclusion can be made that there exists multiple similarity solutions when $\mathrm{S}\ge 2.3224$ and only a single similarity solution exists for$\mathrm{S}<2.3224$ when K = 0.1. However, the range of multiple similarity solution is $2.3769\le \mathrm{S}$ and there also exists only a single solution when $\mathrm{S}<2.3769$ when K = 0.2 (see Figure 4). Furthermore, skin friction increases as suction is increased in the third solution. It is worth mentioning here that there is no range of no solution. This is one of the big reasons that insist us to reconsider and re-examine the whole problem by reproducing all of the results because there are triple solutions in order to provide true knowledge to the readers and researchers.

Figure 5 illustrates the effect of micropolar parameter on the $h’\left(0\right)$. The effect of local couple stress enhanced as the suction increases in the first and third solutions because increasing suction creates additional resistance in the flowing fluid inside the boundary layer. However, increments in the material parameter produce more coupling of stress. Figure 6 shows the nature of the heat transfer rate for various values of the suction. It has been examined that the heat transfer rate increases in the first and second solutions for the higher values of the suction parameter, while the third solution shows opposite compliance.

Finally, we plot Figure 7, Figure 8, Figure 9 and Figure 10 to show the existence of triple solutions of velocity, microrotation, and temperature profiles for different values of material parameter K. In Figure 7, the dual nature of behavior has been noticed in the first solution. Velocity profile increases in the second solution when K is increased; the physical material parameter reduces the effect of drag force due to that thickness of the momentum boundary layer enhanced. On the other hand, the opposite trend has been observed in the third solution. Dual behavior can be noticed microrotation profile in Figure 8 for all solutions. The thickness of thermal boundary layer increases in the first and second solutions as material parameter K is increased, as in Figure 9, since the non-Newtonian parameter produces more viscosity, decreases the velocity of profiles, and forces fluid flow to stay on the hotter surface, as a result temperature of fluid increases and the boundary layer becomes thicker. However, the opposite behavior is noticed for the third solution. The Prandtl effect on the temperature distribution is depicted in Figure 10. It is observed that the temperature of fluid diminishes for the higher values of the Prandtl number for all solutions. Physically, it can be explained, as the Prandtl number ($Pr=\frac{\mu c_p}{k}$) has an inverse relationship with the thermal conductivity and, consequently, diminishes the thickness of the thermal boundary layer.

Table 2. Smallest eigen values $\epsilon$ at several values of S and K when Pr = 1.

K	S	$\mathit{\epsilon }$
		1st solution	2nd solution	3rd solution
0.1	2.3224	1.28061	0	0
-	2.4	1.0662	−0.06382	−0.13406
0.2	2.3769	1.36201	0	0
-	2.4	1.1364	−0.10482	−0.17482

shows the smallest eigenvalues $\epsilon$ for the selected values of S and K. The positive smallest eigenvalue makes the initial disturbance decay and, in this way, the flow becomes stable. Conversely, the negative smallest eigenvalue outcomes in an initial growth of disturbance, in this manner, the flow is unstable. It is seen from Table 2 that $\epsilon$ is negative for the second and the third solutions, while positive for the first solution. Thus, the second and the third solution are not stable, and the first solution is stable. From this discussion, it can be concluded that the first solution of Aurangzaib et al., [29] is not stable and not physically realizable; therefore, in this stage, it could be said that the first solution is actually the second or the third solution.

6. Conclusions

The micropolar fluid over the shrinking surface has been considered. The system of governing equations has been transformed into the system of ODEs by using appropriate exponential similarity transformation. The system of ODEs is reduced to IVPs by employing the shooting method before solving IVPs by the Runge Kutta method. The pointwise conclusions of this study are given below:

• Triple solutions appear.
• There are ranges of multiple solutions and no solutions that depend upon the suction parameter.
• According to stability analysis, the first solution is stable, which can be experimentally seen.
• The results of Aurangzaib et al., [29] are unstable.
• The thickness of thermal boundary layer increases in the first and the second solutions as material parameter K is increased.
• Increments in the material parameter produce more couple stress.