**3. Stability Analysis**

Weidman et al. [42] initiated a study of the stability analysis of multiple solutions. Since then, some researchers, such as Rosca and Pop [43] and Lund et al. [44,45], performed stability analyses in their studies on multiple solutions of fluid flow problems. They found that only the first or upper solution has a physical meaning, while all of the remaining solutions (second or third) are not physically relevant or, in other words, are said to be unstable solutions. The first step in finding the stability of the solutions is to change the momentum, heat, and concentration equations into an unsteady form by considering a new variable τ. In our case, we have τ = *<sup>U</sup>*<sup>0</sup> 2*l e x <sup>l</sup>* ·*t*, as defined in the paper of Rehman et al. [46]:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \left( 9 + \frac{K\_1}{\rho} \right) \frac{\partial^2 u}{\partial y^2} + \frac{\kappa}{\rho} \frac{\partial N}{\partial y} - \frac{9}{K} u - \frac{b}{\sqrt{K}} u^2 - \frac{\sigma B^2 u}{\rho} \tag{18}$$

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

$$\frac{\partial T}{\partial t} + u \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial y} = a \frac{\partial^2 T}{\partial y^2} + \tau\_1 \left[ D\_B \frac{\partial \mathbb{C}}{\partial y} \frac{\partial T}{\partial y} + \frac{D\_T}{T\_{\text{co}}} \left( \frac{\partial T}{\partial y} \right)^2 \right] \tag{20}$$

$$
\rho \frac{\partial \mathbb{C}}{\partial t} + u \frac{\partial \mathbb{C}}{\partial x} + v \frac{\partial \mathbb{C}}{\partial y} = D\_B \frac{\partial^2 \mathbb{C}}{\partial y^2} + \frac{D\_T}{T\_{\text{oo}}} \frac{\partial^2 T}{\partial y^2} \tag{21}
$$

The presence of τ is associated with initial value problems of the stable solution. Equating the derivative of functions with respect to *x* being equal to zero leads to the new similarity transfer variables in the presence of τ and η, which can be expressed as:

$$\begin{array}{lcl}\psi = & \sqrt{2\mathcal{S}l\mathcal{U}l\_0e^{\frac{\pi}{2l}}}f(\eta,\tau);\mathcal{N} = \mathcal{U}e^{\frac{3\pi}{2l}}\sqrt{\frac{\mathcal{U}\_0}{2\mathcal{S}l\mathcal{U}}}g(\eta,\tau);\ \mathcal{O}(\eta,\tau) = \frac{(T-T\_{\text{ov}})}{(T\_w-T\_{\text{ov}})};\\\mathcal{O}(\eta,\tau) = & \frac{(\mathcal{C}-\mathcal{C}\_{\text{ov}})}{(\mathcal{C}\_w-\mathcal{C}\_{\text{ov}})};\ \eta = y\sqrt{\frac{\mathcal{U}\_0}{2\mathcal{N}l}}e^{\frac{\pi}{2l}};\ \tau = \frac{\mathcal{U}\_0}{2l}e^{\frac{\pi}{l}}t\end{array} \tag{22}$$

By applying Equation (22) in Equations (18)–(21), we get:

$$f(1+K)\frac{\partial^3 f}{\partial \eta^3} + f \frac{\partial^2 f}{\partial \eta^2} + K \frac{\partial g}{\partial \eta} - 2 \left(\frac{\partial f}{\partial \eta}\right)^2 - F\_S \left(\frac{\partial f}{\partial \eta}\right)^2 - K\_1 \frac{\partial f}{\partial \eta} - M \frac{\partial f}{\partial \eta} - \frac{\partial^2 f}{\partial \tau \partial \eta} = 0 \tag{23}$$

$$f\left(1+\frac{K}{2}\right)\frac{\partial^2 \mathcal{g}}{\partial \eta^2} + f\frac{\partial \mathcal{g}}{\partial \eta} - 3g\frac{\partial f}{\partial \eta} - 2Kg - K\frac{\partial^2 f}{\partial \eta^2} - \frac{\partial g}{\partial \tau} = 0\tag{24}$$

$$\frac{1}{Pr} \frac{\partial^2 \theta}{\partial \eta^2} + f \frac{\partial \theta}{\partial \eta} + N\_b \frac{\partial \mathcal{Q}}{\partial \eta} \frac{\partial \theta}{\partial \eta} + N\_l \left(\frac{\partial \theta}{\partial \eta}\right)^2 - \frac{\partial \theta}{\partial \tau} = 0 \tag{25}$$

$$\frac{\partial^2 \mathcal{Q}}{\partial \eta^2} + \text{Scf} \frac{\partial \mathcal{Q}}{\partial \eta} + \frac{N\_t}{N\_b} \frac{\partial^2 \theta}{\partial \eta^2} - \text{Sc} \frac{\partial \mathcal{Q}}{\partial \tau} = 0 \tag{26}$$

subject to the boundary conditions:

$$f(0, \tau) = \text{S; } \frac{\partial f(0, \tau)}{\partial \eta} = -1 + \lambda \frac{\partial^2 f(0, \tau)}{\partial \eta^2}; \text{g}(0, \tau) = -m \frac{\partial^2 f(0, \tau)}{\partial \eta^2}; \theta(0, \tau) = 1; \mathcal{Q}(0, \tau) = 1 \tag{27}$$

$$\frac{\partial f(\eta, \tau)}{\partial \eta} \to 0; \text{ g}(\eta, \tau) \to 0; \ \theta(\eta, \tau) \to 0; \mathcal{Q}(\eta, \tau) \to 0 \quad \text{as } \eta \to \infty$$

In order to indicate the solution stability of *f*(η) = *f*0(η), *g*(η) = *g*0(η), θ(η) = θ0(η) and ∅(η) = ∅0(η), which satisfy the equation of the boundary value problem (23)–(26) with boundary condition (27), we follow the suggestion of Rehman et al. [46] by introducing the following functions:

$$\begin{cases} f(\eta,\tau) = f\_0(\eta) + e^{-\varepsilon\tau} F(\eta,\tau) \\ g(\eta,\tau) = g\_0(\eta) + e^{-\varepsilon\tau} G(\eta,\tau) \\ \theta(\eta,\tau) = \theta\_0(\eta) + e^{-\varepsilon\tau} H(\eta,\tau) \\ \mathcal{Q}(\eta,\tau) = \mathcal{Q}\_0(\eta) + e^{-\varepsilon\tau} S(\eta,\tau) \end{cases} \tag{28}$$

where *F*(η, τ), *G*(η, τ), *H*(η, τ) and *S* (η, τ) are small relative to *f*0(η), *g*0(η), θ0(η) and ∅0(η) of the steady state solutions. It should be noted that the range of these functions are 0 < *F*(η, τ) < 1, 0 < *G*(η, τ) < 1, 0 < *H*(η, τ) < 1 and 0 < *S*(η, τ) < 1. Furthermore, ε is an unknown eigenvalue parameter, which needs to be found. Substituting the values of the functions and their derivatives from Equation (28) in Equations (23)–(26) with the boundary condition (27), we have:

$$(1+K)F\_{0}^{\prime\prime} + f\_{0}F\_{0}^{\prime\prime} + F\_{0}f\_{0}^{\prime\prime} + KG\_{0}^{\prime} - 4f\_{0}^{\prime}F\_{0}^{\prime} - 2F\_{S}f\_{0}^{\prime}F\_{0}^{\prime} - K\_{1}F\_{0}^{\prime} - MF\_{0}^{\prime} + \varepsilon F\_{0}^{\prime} = 0 \tag{29}$$

$$\left(1+\frac{K}{2}\right)G\_0'' + f\mathfrak{g}G\_0' + F\mathfrak{g}G\_0' - 3\mathfrak{g}\mathfrak{a}F\_0' - 3\mathfrak{g}\mathfrak{a}F\_0' - 2K\delta\mathfrak{G}\_0 - K\delta F\_0''' + \varepsilon\mathfrak{G}\_0 = 0\tag{30}$$

$$\frac{1}{Pr}H\_0'' + f\_0H\_0' + F\_0\theta\_0' + Nb\mathcal{Z}\_0'H\_0' + NbS\_0'\theta\_0' + 2Nt\theta\_0'H\_0' + \varepsilon H\_0 = 0\tag{31}$$

$$\mathcal{S}\_0^{\prime\prime} + \text{Sc}\Big(f\_0 \mathcal{Q}\_0^{\prime} + F\_0 \mathcal{S}\_0^{\prime}\Big) + \frac{\text{Nt}}{\text{Nb}} H\_0^{\prime\prime} + \text{Sc} \cdot \varepsilon \mathcal{S}\_0 = 0\tag{32}$$

subject to the boundary conditions:

$$\begin{aligned} F\_0(0) &= 0, F\_0'(0) = \lambda F\_0''(0), G\_0(0) = -m F\_0''(0), \; H\_0(0) = 0, \; S\_0(0) = 0 \\ F\_0'(\eta) &\to 0, \; G\_0(\eta) \to 0, H\_0(\eta) \to 0, S\_0(\eta) \to 0, \; \text{as } \eta \to \infty \end{aligned} \tag{33}$$

We assumed τ = 0 for Equations (23)–(26) in order to calculate the initial growth and decay of the solution of Equation (28), as recommended by Alarifi et al. [47]. Under these circumstances, *F*(η, τ), *G*(η, τ), *H*(η, τ) and *S* (η, τ) can be written as *F*0(η), *G*0(η), *H*0(η) and *S*0(η), respectively.

It is stated in the studies of Lund et al. [44,45] and Haris et al. [48] that eigenvalues can be determined if and only if the boundary condition of any one function of the following functions *F*0(η), *G*0(η), *H*0(η) and *S*0(η) can be relaxed into the initial condition by differentiating that function one more time. In this study, we relaxed *F*0(η) → 0 as η → ∞ and then solved the system of Equations (29)–(32) subject to Equation (33) along with the new relaxed boundary condition *F*-- <sup>0</sup>(0) = 1. It is worth mentioning that the sign of the smallest eigenvalues (ε) determines the stability of the solutions. The smallest eigenvalue is negative (positive), which indicates that the solution of the flow is unstable (stable) and that there is an initial growth (decay) of disturbances.

#### **4. Results and Discussion**

In order to fully understand the considered fluid flow model, the numerical study has been carried out for various important physical parameters, such as the magnetic parameter *M*, permeability parameter *K*1, Forchheimmer parameter *FS*, non-Newtonian parameter *K*, thermophoresis parameter *Nt*, Brownian motion parameter *Nb*, etc. The highly non-linear system of the quasi-ordinary differential Equations (12)–(15), along with the boundary conditions (16), have been solved by using the shooting

method, and triple solutions were found. The value of η<sup>∞</sup> is chosen from 4 to 8, and it is worth noting that the value of η<sup>∞</sup> increases until the profiles of the velocity, temperature, and concentration converge asymptotically to the momentum, temperature, and concentration boundary layers, respectively.

Figure 2 was drawn to analyze the effect of the material parameter *K* on the velocity profiles. The thickness of the velocity boundary layer increases when the micropolar parameter *K* is increased in the first and the second solutions, due to the fact that the material parameter reduced the drag force and that the hydrodynamic boundary layer therefore rose. We noted that when *K* = 0 (Newtonian fluid), the second solution did not exist. On the other hand, the dual nature of the velocity profile was observed in the third solution. Figure 3 depicts the variation of the velocity profiles for the different values of the Forchheimmer parameter *FS*. We observed that, due to increments in *FS*, the resistant force occurred when the fluid flow was flowing on the porous surface, and hence the velocity of the flow declined in the third solutions. The dual behavior of the velocity profile was noticed in the second solution. However, no change could be seen in the first solution when *FS* was increased. The velocity and thickness of the momentum boundary layers are inversely (directly) proportional to *K*<sup>1</sup> and *M* in the first (third) solution. On the other hand, the dual behavior of the velocity profile was seen in the second solution, as illustrated in Figures 4 and 5. The effect of the slip parameters λ and *m* on the velocity profiles are shown in Figures 6 and 7. In the first solution, the velocity boundary thickness was reduced as λ and *m* increased; this was due to the fact that the velocity of the fluid and surface have a big difference when the velocity slip factor is enhanced. For the second (third) solution, the velocity profiles decreased (increased) initially after they inclined (declined) when λ and *m* improved. The dual nature of the flow in some sense indicates that there is an initial growth of disturbance. Initially, the microrotation profile was reduced and then started to rise with increasing values of the micropolar parameter in the first solution, as shown in Figure 8. The thickness of the microrotation boundary layer was enhanced (reduced) as *K* was enhanced in the third (second) solution. The dual behavior of the microrotation profile was been observed in all solutions except the first solution when the value of *m* increased. The thickness of the microrotation boundary layer increased with increasing values of *m* in the first solution, as demonstrated in Figure 9. Figure 10 was drawn to analyze the variation of the temperature profiles for different values of the Prandtl number *Pr*. The thermal boundary layer thickness and temperature were incrementally reduced in the values of the Prandtl number *Pr* for all of the solutions, as expected. This is due to the fact that a high Prandtl number causes the thermal conductivity of nanofluid to diminish, and as a result the temperature is reduced. The variation of the temperature profiles for different values of the Brownian motion parameter *Nb* is shown in Figure 11. It was observed in all solutions of the temperature profiles that the temperature and thermal layer thickness were enhanced with increasing values of *Nb*. This is due to fact that the Brownian motion *Nb* increases the kinetic energy of the nanofluid; thus, the temperature of the nanofluid increases. The temperature profile, with an effect resulting from the thermophoresis parameter *Nt*, is illustrated in Figure 12. In all three solutions of the nanofluid flow problem, we noted that as the thermophoresis parameter *Nt* rose, the temperature and thickness of the thermal layer increased. This is because the thermophoretic force is generated by *Nt* and the temperature gradient, which pushes the flow of the nanofluid far from the boundary layer as a resulting thickness of the thermal boundary layer increases. The thickness of the concentration boundary layer declined with increasing values of the Brownian motion parameter *Nb* in all solutions in Figure 13, which was expected. This was physically justified by the fact that Brownian motion is generated when nanoparticle and base fluid are mixed together in a nanofluid system. Since Brownian diffusion shows the conduction of heat under those circumstances, the thickness of the concentration boundary layer decreases. Figure 14 was sketched to examine the effect of *Nt* on the concentration profile of nanoparticles. In all three solutions, the thickness of the concentration boundary layer was enhanced when thermophoresis increased.

**Figure 2.** Profiles of the velocity for increasing values of *K*.

**Figure 3.** Profiles of the velocity for increasing values of *FS*.

**Figure 4.** Profiles of the velocity for increasing values of *K*1.

**Figure 5.** Profiles of the velocity for increasing values of *M*.

**Figure 6.** Profiles of the velocity for increasing values of λ.

**Figure 7.** Profiles of the velocity for increasing values of *m*.

**Figure 8.** Microrotation profiles for increasing values of *K*.

**Figure 9.** Microrotation profiles for increasing values of *m*.

**Figure 10.** Temperature profiles for increasing values of *Pr*.

**Figure 11.** Temperature profiles for increasing values of *Nb*.

**Figure 12.** Temperature profiles for increasing values of *Nt*.

**Figure 13.** Nanoparticle concentration profiles for increasing values of *Nb*.

**Figure 14.** Nanoparticle concentration profiles for increasing values of *Nt*.

Wang [49] and many researchers also stated in their studies that the similarity solution of the fluid flow problems over the shrinking surface is possible to obtain when sufficient wall mass suction is applied. The flow of Newtonian fluid is different from that of non-Newtonian fluid; it is observed for micropolar nanofluid that when the value of the micropolar parameter *K* increases, a strong mass suction is required to obtain the solution. In this study, we discovered that there exist two regions for similarity solutions, namely multiple solutions and single solution, depending on the mass suction parameter. For *K* = 0.1(*K* = 0.2) there is a range of triple solutions when *S* ≥ 2.0337(*S* ≥ 2.7148), and a single similarity solution exists, *S* < 2.0337(*S* < 2.7148), as shown in Figure 15. When the suction is enhanced, the skin friction increases in the first solution and decreases in the second and the third solutions. Figures 16–18 were drawn to examine the effect of the suction *S* and micropolar parameter *K* on the couple stress coefficient, and the heat and concentration transfer rate, respectively. In all graphs, when the suction is increased, the couple stress coefficient, heat, and concentration transfer rate are enhanced for all of the solutions.

**Figure 15.** Coefficient of the skin friction with various values of *S*.

**Figure 16.** Couple stress coefficient with various values of *S*.

**Figure 17.** Heat transfer rate with various values of *S*.

**Figure 18.** Concentration transfer rate with various values of *S*.

By performing a stability analysis, the stability of the fluid flow solutions is achieved in this research. We need to perform a stability analysis when more than one solution exists in the flow problem. The main focus of this analysis is to determine which solution is stable and physically possible and which one is unstable. It is worth noting that the stability of the solution depends on the sign of the smallest eigenvalue. The value of the smallest eigenvalue is determined through Equation (24), for which we have to solve Equations (25)–(28), along with the boundary conditions (29). The smallest eigenvalues ε are demonstrated in Table 1 for different values of the suction and non-Newtonian parameters. It is clear from Table 1 that the negative (positive) values of the smallest eigenvalue ε

indicate an initial growth (decay) of the disturbance that will interrupt (resume) the boundary layer separation and flow from becoming unstable (stable). It is worth mentioning that the stable solution always provides a good physical meaning which can be realized.


**Table 1.** Smallest eigenvalues for different values of *K* and *S* when *M* = 0.5, *Pr* = *Sc* = 1, *m* = 0.5, *K*<sup>1</sup> = 0.1, *FS* = 0.2, λ = 0.1, *Nb* = 0.2 and *Nt* = 0.15.
