**6. Results and Discussion**

The mathematical model in Equations (23)–(25) was solved numerically by means of the boundary value problem solver function bvp4c in the MATLAB software. The numerical results were derived while limiting the relative tolerance to 1 × 10−10. Some of the governing parameter values were fixed throughout the computation process to align with the motivation of this study. For example, the Prandtl

number (Pr) value was fixed at 4.36 because it represents the base fluid, *n*-hexane. The power-law index (*<sup>m</sup>*0) was fixed at 1.5 to investigate the dilatant features of the Sutterby fluid. The obtained non-uniqueness solutions were classified based on how early the solution converged asymptotically. For example, the numerical solution that converged earlier asymptotically in the velocity/temperature profiles was labelled the first solution. The other solution that converged later, asymptotically, was labelled as the second solution. Before presenting the numerical results, we provide validation of our numerical method by solving the model presented in [41] and compare the numerical results with the results reported by [41]. Table 2 shows the comparison results, and there is a good agreement. Bhattacharyya et al. [41] employed the shooting method to solve the model, and Table 2 confirms that the bvp4c function is capable of precisely solving the boundary value problem.


**Table 2.** Numerical validation of *f* (0) when *S* = 0 in [41].

\**c*/*a* is the stretching/shrinking parameter in [41].

Figure 2 shows the influence of the suction parameter (*s*) on the reduced skin friction coefficient -*Cf x*Re1/2 *x* and velocity profiles (*f*(η)). Based on the first solution in Figure 2a, an increment in *s* increases the values of *Cf x*Re1/2 *x* past a shrinking sheet. Primarily, an increment in *s* from 6 to 9 strengthens the impact of suction at the surface of the shrinking sheet. The act of suction encourages the laminar flow by trapping the low speed fluid molecules in the boundary layer region. This then leads to increasing of the fluid velocity past the shrinking sheet and is illuminated in Figure 2b. The increment of the fluid velocity reduces the momentum boundary layer thickness and increases the wall shear stress over the shrinking sheet. The high wall shear stress eventually increases the values of *Cf x*Re1/2 *x* as *s* increases. Interestingly, the second solution in Figure 2a shows the opposite trend to the first solution, where an increment in *s* decreases the values of *Cf x*Re1/2 *x* . The state of suction, which was interpreted as enhancing the fluid velocity, is now seen to decrease the fluid velocity and increase the momentum boundary layer thickness (see Figure 2). The saturated state of the shrinking sheet may be the cause of these consequences. Later, the reducing fluid velocity lowers the wall shear stress and then decreases the values of *Cf x*Re1/2 *x*as *s* increases.

Figure 3 demonstrates the impact of the Deborah number (*De*) on *Cf x*Re1/2 *x* and velocity profiles. Both solutions in Figure 3a convey that an increment in *De* augments the values of *Cf x*Re1/2 *x* as the sheet is shrinking. The Deborah number is used to enlighten the viscoelastic feature of a material. Here, an increment in *De* results in the increment of the shear thickening Sutterby fluid velocity. The valuable work of Azhar et al. [34] also reported a similar trend. The increment of the fluid velocity then enhances the wall shear stress past the shrinking sheet and increases the values of *Cf x*Re1/2 *x* . The increment in *s* and *De* assist in delaying the flow separation significantly.

**Figure 2.** Impact of the suction parameter (*s*) on: (**a**) the reduced skin friction coefficient; (**b**) velocity profiles as *s* varies when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36, *De* = 1.5, *M* = 0.5, and φ = 0.02.

**Figure 3.** Impact of the Deborah number (*De*) on: (**a**) the reduced skin friction coefficient; (**b**) velocity profiles as *De* varies when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36,*s* = 7, *M* = 0.5, and φ = 0.02.

The first and second solutions in Figure 4a lead to a decrement in *Cf x*Re1/2 *x* when *M* increases from 0.5 to 1.0. The magnetic field presents in an electrically conducting fluid as an electromagnetic force in the fluid flow region, which slows the fluid moving past the shrinking sheet. This is reflected by the velocity profiles in Figure 4b, where the fluid velocity decreases when *M* increases. The decrement in the fluid velocity then leads to an increase of the momentum boundary layer thickness and decreases the wall shear stress past the permeable shrinking sheet. Thus, the values of *Cf x*Re1/2 *x* decrease with the rising value of *M*. Unlike the shrinking case, different fluid flow behavior is perceived in the first solution when the Sutterby nanofluid flows towards a permeable stretching sheet. When the magnetic effect increases past a stretching sheet, the fluid velocity increases, although the increment is not significant. This is agreeable because the fluid flow is in the same direction as the stretching sheet and the action of the stretching sheet speeds up the fluid flow. The increment in fluid velocity reduces the momentum boundary layer thickness, increases the wall shear stress, and enhances the value of *Cf x*Re1/2 *x* . The negative values of *Cf x*Re1/2 *x* indicate that the stretching sheet imposes a drag force on the fluid. Moreover, the reverse flow is observed through the second solution's presence when the permeable sheet is stretching in Figure 4a. The velocity profiles in Figure 4c support this by displaying the velocity overshoot (see the second solution profiles) when *M* varies. Thus, it is clear that reverse

flow does exist in the stretching sheet case, and this may be due to the state of the sheet where the suction intensity is weak when the effect of *M* increases.

**Figure 4.** Impact of the magnetic parameter (*M*) on: (**a**) the reduced skin friction coefficient; (**b**) velocity profiles as *M* varies past a shrinking sheet (ε = −<sup>4</sup>); (**c**) velocity profiles as *M* varies past a stretching sheet (ε = 7) when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36, *De* = 1.5,*s* = 7, and φ = 0.02.

Figure 5 portrays the effect of the nanoparticle volume fraction or φ on *Cf x*Re1/2 *x* and velocity profiles. The increment in φ increases the values of *Cf x*Re1/2 *x* over a permeable shrinking sheet. An increased ratio of φ in the base fluid increases fluid viscosity, which then enhances the fluid velocity past the permeable shrinking sheet (see Figure 5b). These then affect the wall shear stress to increase and, consequently, raise the values of *Cf x*Re1/2 *x* . Velocity overshoots in the boundary layer are apparent in Figures 2b, 3b, 4b and 5b. These velocity overshoots near the permeable shrinking sheet indicate that the fluid velocity is higher than the shrinking sheet's velocity [42].

Table 3 exhibits the effect of the radiation parameter (*Rd*) on the reduced local Nusselt number -Re−1/2 *x Nux x*<sup>−</sup>2/5 over the permeable shrinking surface. Both solutions allude to the enhancement of Re−1/2 *x Nux x*<sup>−</sup>2/5 when the impact of radiation grows in the fluid flow region. An increment in *Rd* hints at the release of energy in the form of heat from the fluid flow and decreases the fluid temperature profile. Thus, the thermal boundary layer becomes thinner and the wall heat flux increases. The depreciation in the thermal conductivity induces an increase in the rate of heat transfer or -Re−1/2 *xNux <sup>x</sup>*<sup>−</sup>2/5.

Furthermore, the magnetic parameter and the nanoparticle volume fraction have minimal effect in delaying flow separation. This is evident by the critical values (<sup>ε</sup>0), as shown in Figures 4a and 5a.

**Figure 5.** Impact of the nanoparticle volume fraction (φ) on: (**a**) the reduced skin friction coefficient; (**b**) velocity profiles as φ varies when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36, *De* = 1.5,*s* = 7, *M* = 0.5, and φ = 0.02.

**Table 3.** The effect of the radiation parameter (*Rd*) on the reduced local Nusselt number when *M* = 0.5, *m*0 = 1.5, Pr = 4.36, φ = 0.02,*s* = 7, and *De* = 1.5 as ε varies.


The results of the stability analysis are presented in Table 4. The first solution achieves the positive eigenvalues while the second solution attains the negative eigenvalues. Based on the signs of eigenvalues, one can say that the positive eigenvalues specify the first solution as a stable solution; the stable solution can be understood as feasible and able to overcome the growth of an initially given disturbance. Furthermore, the negative eigenvalues reveal the second solution as an unstable solution associated with flow separation. The second solution promotes the growth of an initially given disturbance and hence achieves the negative eigenvalue. However, it is vital to identify and verify the stability of non-unique solutions so that the variety of possibilities of fluid flow behavior can be predicted.

Figure 6 depicts the streamlines of the Sutterby fluid under a number of settings. In particular, Figure 6a shows the streamlines when the sheet is impermeable and stretching at the rate of 1.4, while Figure 6b illustrates the streamlines when the sheet is impermeable and shrinking. The reverse flow in Figure 6b is noticeable and proves that the shrinking sheet's state instigates the reverse flow. Next, the streamlines for the fluid flow under the suction influence can be examined (see Figure 6c,d). The reverse flow is now absent past the permeable shrinking sheet (Figure 6d). Thus, it is proved that mass suction succeeds in sustaining the laminar boundary layer flow over a shrinking surface. Figure 6e,f shows the behavior of fluid flow when the rate of stretching or shrinking increases; the fluid pattern being pulled at the surface of the sheet is clear, and again the reverse flow is absent.

**Table 4.** Lowest eigenvalues (γ1) when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36, *De* = 1.5,*s* = 7, *M* = 0.5, and φ = 0.02 as ε varies.

**Figure 6.** Streamlines when *Rd* = 1.2, *m*0 = 1.5, Pr = 4.36, *De* = 1.5, φ = 0.02, *M* = 0.5; (**a**)*s* = 0, ε = 1.4; (**b**) *s* = 0, ε = −1.4; (**c**) *s* = 7, ε = 1.4; (**d**) *s* = 7, ε = −1.4; (**e**) *s* = 7, ε = 4; (**f**) *s* = 7, ε = −4.
