**4. Results and Discussion**

The mathematical computations of this research work were achieved by employing the bvp4c function in the MATLAB programming system subject to the governing ordinary di fferential Equations (7) and (8) together with the boundary conditions (9). A bvp4c method is a notable tool for solving the boundary value problem that has been extensively established by various researchers to clarify the boundary value concern. In order to ensure the obtainment of desired solutions, early estimation of the primary mesh point and variations step size is crucial. Also, a reasonable assumption of the thickness of the boundary layer along with e ffective preliminary approximation is essential for defining the non-uniqueness solutions. The consistency of the results generated in the present study is evaluated with those in [34,58], as accessible in Table 3, which is in excellent agreement.

In this study, we recognized the hybrid Al2O3-Cu/H2O nanofluid by dispersing the first nanoparticle φ1 (alumina) into the base fluid (water) followed by the second nanoparticle φ2 (copper) with various amounts of volume fractions. It is important to declare that the dispersion of discrete nanoparticle Al2O3/Cu is capable of developing Al2O3-Cu/H2O and Cu-H2O nanofluids. Apart from that, the Prandtl (Pr) number is set to be fixed at Pr = 6.2 corresponds to water as the reference-based fluid, while as for the hybrid Al2O3-Cu/H2O nanofluid, the size of the nanoparticles is assumed to be standardized, and the thermophysical properties e ffect of nanoparticles agglomeration is ignored. The main purpose of the present study is to examine the influence of the control parameter such as the nanoparticles volume fraction (φ1, φ2), the unsteadiness parameter (ε), the velocity slip parameter (γ), and the Biot number (Bi) towards the coefficients of skin friction variations (*f* (0)) and the heat transfer rate (θ(0)). The alumina and copper nanoparticles volume fraction in the present work is chosen within the range of 0.00 ≤ φ2 ≤ 0.04 which motivated by the experimental work done by Suresh et al. [65] who conducted the synthesis, characterization of Al2O3 − Cu/H2O nanocomposite powder for different volume concentrations 0.1%, 0.33%, 0.75%, 1%, and 2%. In their valuable study, the stability of the prepared nanofluids was determined by measuring the pH of nanofluids, and the nanofluid stability was found to diminish with increasing volume concentration. On another note, various values of the controlling parameter that has been used, are set within the following extent; 0.1 ≤ ε ≤ 0.2, 0.1 ≤ γ ≤ 0.4, and 0.2 ≤ Bi ≤ 0.7 to ensure the certainty of the obtained solutions. The present study also interested in witnessing the non-uniqueness solutions that appear in the ordinary governing differential of the ensuing problem, hence confirm the real and valid solutions through the stability analysis by utilizing the bvp4c approach in MATLAB operating system (MATLAB R2019b, MathWorks, Natick, MA, USA).



The generated results of Equations (7) and (8) perceive the non-uniqueness (dual) solutions together with the boundary conditions (9) to a specific scope of λ*c* where λ*c* manifests the non-uniqueness solutions meeting point, and this too is regarded as a critical point. There is a single solution at a critical point, and that is the distinct line. The flow separation happens to occur after the critical point, and the flow is no longer laminar and does not obey the principle of boundary layer theory. According to the numerical results attained in this research work, it is proven that a non-uniqueness solutions appear, namely first and second solutions, as depicted in Figures 2–15. The existence of the non-uniqueness solutions contributes to the analysis of solution stability so that one may be able to verify the theoretically relevant solution. The first solution, however, is predicted to be reliable and fundamentally exist in practice. The smallest eigenvalues ω1 for some values of λ when ε = 0.1,φ1 = φ2 = 0.02, Bi = 0.2, and Pr = 6.2 with γ = 0.1, 0.2, 0.4, are tabulated in Table 4. The values of ω1 is noted approaching zero in the first and second solutions for λ → λ*c* in selected cases of γ. Hence, the authors may assure that the stability formulation and practices of the current problem are accurate and reliable. The flow is measured unstable if ω1 is negative because it implies an early development of disruptions that leads to flow separation. The smallest eigenvalue ω1 with positive value implies the flow is practically attainable and sustainable, which elucidates the solution stabilizing property to overwhelm the allowing disruptions. Also, it denotes an initial deterioration of disruptions that appear. The positive value of ω1 describes the stable mode of the flow as in the first solution. The significance of dual solutions that contribute to stability analysis is in revealing the other possibilities of flow behavior, which may be useful for future references in the extrusion process. For example, in the present study, we have presented that the second solutions are unstable solutions with negative eigenvalues. These unstable solutions with the negative eigenvalues infer the growth of the disturbance in the solutions, especially at some range of the stretching/shrinking parameter (λ < <sup>−</sup>1.0). When the shrinking rate increases, it limits other external forces' effect at the sheet and

eventually showed the opposite behavior of the transport phenomena than the first solutions in the same range of λ < −1.0, as the respective parameter varies. Moreover, any boundary value problem can generate more than one solution because of the nonlinearity in the boundary value problem; for example, see the mathematical model (7)–(9). Also, changes in the governing parameter values cause bifurcations in solutions that yield the existence of dual solutions [1]. Therefore, the mathematical model (7)–(9), which obeyed the boundary layer assumptions, managed to exhibit the variety in the fluid flow and heat transfer behavior through the uniqueness and existence of the solutions.

**Figure 2.** Variants of *f* (0) towards λ with γ = 0.1, 0.2, 0.4.

**Figure 3.** Variants of −θ(0) towards λ with γ = 0.1, 0.2, 0.4.

**Figure 4.** Velocity profiles of *f*(η) with λ = −1.3 (shrinking case).

**Figure 5.** Temperature profiles of <sup>θ</sup>(η) with λ = −1.3 (shrinking case).

**Figure 6.** Variants of *f* (0) towards λ with φ2 = 0.00, 0.02, 0.04.

**Figure 7.** Variants of −θ(0) towards λ with φ2 = 0.00, 0.02, 0.04.

**Figure 8.** Velocity profiles of *f*(η) with λ = −1.4 (shrinking case).

**Figure 9.** Temperature profiles of <sup>θ</sup>(η) with λ = −1.4 (shrinking case).

**Figure 10.** Variants of *f* (0) towards λ with ε = 0.10, 0.15, 0.20.

**Figure 11.** Variants of −θ(0) towards λ with ε = 0.10, 0.15, 0.20.

**Figure 12.** Velocity profiles of *f*(η) with ε = 0.10, 0.15, 0.20.

**Figure 13.** Temperature profiles of <sup>θ</sup>(η) with ε = 0.10, 0.15, 0.20.

**Figure 14.** Variants of −θ(0) towards λ with Bi = 0.2, 0.5, 0.7.

**Figure 15.** Temperature profiles of <sup>θ</sup>(η) with Bi = 0.2, 0.5, 0.7.


**Table 4.** Smallest eigenvalues ω1 for some values of λ when γ = 0.1, 0.2, 0.4.

The influence of velocity slip on the skin friction coe fficient and the local Nusselt number past a convectively heated stretching/shrinking sheet of hybrid Al2O3-Cu/H2O nanofluid are displayed in Figures 2–5. Figure 2 presents the coe fficient of skin friction (*f* (0)) towards the stretching/shrinking parameter λ, which revealed a decline of *f* (0) with the occurrence of velocity slip e ffect (γ = 0.1, 0.2, 0.4) at the boundary. The same results are obtained in the previous literature, as reported by Dzulkifli et al. [59]. Based on a physical perspective, an improvement in the slip parameter reflects the fact that the vorticity produced by the stretching/shrinking velocity is gradually decreased, thus the vorticity stays restricted within the boundary layer for greater stretching/shrinking velocity with the same straining velocity of the stagnation flow, and subsequently, the steady solution is achievable for some broad values of λ, as stated by Mahapatra and Nandy [34]. Figure 2 emphasizes that the solution for a certain value of γ persist prior to a critical value λ = (<sup>λ</sup>*c* < 0) at which the boundary layer splits from the convectively heated stretching/shrinking sheet and the solution on the basis of the boundary layer approximations is not feasible. In short, no solution exists for λ*c* < −1.6851. Moreover, the expansion of γ results in the increment of |λ*c*| suggesting that the velocity slip parameter is e fficiently influenced by raising the range of dual solutions. Thus, it is proven that the existence of a slip velocity impact may prolong the separation of the boundary layer. Figure 2 also highlights as the sheet is stretching at the rate of λ = 1, the value of *f* (0) = 0 which describes no frictional drag exerted at the convectively heated stretching/shrinking sheet.

Figure 3 exposes an upward trend in −θ (0) when the velocity slip parameter arises on the convectively heated stretching/shrinking sheet, which is proportionate to the heat transfer rate. Those findings are in line with the remarkable work reviewed by Mahapatra and Nandy [34] and Khashi'ie et al. [66]. Evidently, in the case of the first solution, the rate of heat transfer presents an ascendant trend with an increase of γ while the second solution permits a top-down direction when the velocity slip γ enlarges. From the existing and current findings, the authors can conclude that the velocity slip contributes to the improvement of the heat transfer rate significantly, prior to this case study. However, the authors would also like to declare that such results may vary if di fferent control parameters are taken into consideration. The distribution of velocity (*f* (η)) and temperature (θ(η)) profiles over several values of γ in the case of a convectively heated shrinking sheet is certified in Figures 4 and 5. As demonstrated in Figure 4, it is confirmed that as γ increases, *f* (η) is subsequently decreased in the first solution and revealed an upward trend in the velocity profile in the second solution with the presence of velocity slip condition. As the slip occurs, the velocity flow near the sheet in no longer equal to stretching velocity. This is due to the slip condition where the pulling force of the stretching sheet is partly shifted to the fluid, resulting in the decrement of the velocity profiles in the first solution. In short, the slip condition reduces the momentum transfer from the sheet to the fluid. A contradict results are obtained in <sup>θ</sup>(η) as the velocity slip e ffect is exaggerated, as shown in Figure 5. The first solution in the temperature profile distribution increases as γ intensifies, while it diminishes in the second solution. This may occur due to the appearance of slip boundary conditions on the wall, which triggers the hybrid nanofluid to retain the velocity on the walls, and such slip may prevent from total heat exchange of the hybrid nanofluid. It is therefore noticed when the slip coefficient is applied, the difference in temperature rate is increased. Additionally, the effect of convective heating progress that has been reflected in the present study, which controls the temperature of the stretching/shrinking surface also might contribute to this phenomenon. As a matter of fact, greater convection leads to increased surface temperatures, which permit the thermal impact to penetrate deeper within the hybrid Al2O3-Cu/H2O nanofluid.

Figures 6 and 7 expose the coefficient of skin friction (*f* (0)) and the rate of heat transfer (−θ(0)) of a conventional nanofluid (φ1 = 0.02, φ2 = 0) and hybrid nanofluid (φ1 = 0.02,φ2 = 0.02, 0.04) past a convectively heated stretching/shrinking sheet when φ2 varies from 0.00 to 0.02 where ε = 0.1, γ = 0.2, Bi = 0.2, and Pr = 6.2. Figure 6 manifests that addition in φ2 which indicates the transformation of the conventional Al2O3-H2O fluid to the hybrid Al2O3-Cu/H2O nanofluid, upsurges the values of *f* (0) once the sheet is shrinking. The viscosity of hybrid Al2O3-Cu/H2O nanofluid rises when φ2 increases, which eventually improves the fluid velocity over the convectively heated shrinking sheet, as proven in Figure 8. The velocity profile in Figure 8 clarifies that the momentum boundary layer thickness was diminished in response to the rise of φ2, thereby raising the velocity of the fluid and boosting the gradient of velocity. In fact, the thinner momentum boundary layer continues to evolve the wall shear stress as well as the convectively heated shrinking sheet, leading to an enhancement of *f* (0). As *f* (0) increases, the result implies the increase of the frictional drag exerted on the convectively heated shrinking surface, which may delay the boundary layer flow separation. Besides, Figure 6 also highlights when the sheet is stretching at the rate of λ = 1, the value of *f* (0) = 0 which explains no frictional drag exerted at the sheet surface. Meanwhile, Figure 7 illustrates an increasing trend of the heat transfer characteristic or −θ(0) when the values of φ2 past a convectively heated shrinking sheet, and this trend holds true to the first solution but is in dispute with the results in the second solution. In essence, as the conventional Al2O3-H2O nanofluid becomes the hybrid Al2O3-Cu/H2O nanofluid, the heat transfer efficiency improves. This finding upholds the assumption of the convective heat transfer system can be improved by optimizing the nanoparticle concentration when φ2 increased. The results obtained in Figures 6 and 7 are consistent with Waini et al. [20] and Zainal et al. [21], whereby adding the concentrations of hybrid nanoparticles may contribute to the improvement of the heat transfer rate, accordingly. The temperature profile in Figure 9 describes the temperature variations when the conventional Al2O3-H2O nanofluid becomes the hybrid Al2O3-Cu/H2O nanofluid in both first and second solutions. The incline in the temperature of hybrid nanofluid proliferates the thermal conductivity, which may be triggered by the extra energy dispersed through the increment of the nanoparticles volume fraction over the state of convectively heated stretching/shrinking sheet.

Figures 10–13 show the impact of the unsteadiness parameter (ε) towards a convectively heated stretching/shrinking sheet when ε shifts from 0.1 to 0.2. The hybrid Al2O3-Cu/H2O nanofluid characteristic with regard to the coefficient of skin friction (*f* (0)) is depicted as in Figure 10. Figure 10 captures that when the sheet shrinks, the increment in ε conclusively increases the trend of *f* (0) in the first solution. An increment in the unsteadiness parameter results in the reduction of the boundary layer thickness, as depicted in Figure 12 and consequently upsurge the velocity gradient on the convectively heated stretching/shrinking sheet, thus *f* (0) improves. Furthermore, the existence of nanoparticle volume fraction in the working fluid (Al2O3-Cu/H2O) might also trigger the increment of *f* (0) owing to an uplift of the hybrid nanofluid viscosity. This result is aligned with the study done by Ismail et al. [67]. Figure 11 presents the unsteadiness parameter effect towards the rate of heat transfer <sup>−</sup>θ(0). According to the generated results, the heat transfer rate increases when ε increases when the convectively heated sheet is shrinking. The forced convective heat transfer is indeed proportional to the effectiveness of nanofluid heat conductivity, hence raising the reduced local Nusselt number notably. The dimensionless velocity profiles *f*(η) with a different value of ε are depicted in Figure 12, where the presence of dual velocity profiles is also observed. As illustrated in Figure 12, the first solution increases proportionally to the increment of ε values while the second

solution displayed contradict results of the first solution, possibly because of the enhancement in the unsteadiness of the flow. Meanwhile, the same trend of the solution in Figure 12 also reflected the graph of temperature distribution <sup>θ</sup>(η) in the convectively heated stretching/shrinking sheet with the existence of unsteadiness parameter, as portrayed in Figure 13. Apart from that, it is proven that the second solution in both profiles, i.e., velocity and temperature distribution showed larger boundary layer thickness in each solution of the unsteady cases than those of the first solution.

Figures 14 and 15 depict the variants of heat transfer rate −θ(0) and temperature distribution profile <sup>θ</sup>(η) with a different value of the Biot number (Bi = 0.2, 0.5, 0.7) towards the convectively heated stretching/shrinking sheet. In Figure 14, the heat transfer rate shows an upsurge trend in the dual solutions along with the augmentation of Bi values. The Biot number signifies the conduction resistance ratio within the sheet to convection resistance at the sheet. The critical values of the different usage of Biot numbers sugges<sup>t</sup> no substantial impact on the magnitude of <sup>−</sup>θ(0), as highlighted by Jusoh et al. [68]. An increment in the Biot number related to the improvement of convective heating is observed to decrease the fluid temperature proficiently in the first and second solutions, as displayed in Figure 15. This result is in contrast with the idea of a large Biot number representing larger internal thermal resistance of the sheet compare to the boundary layer thermal resistance since the temperature distribution profile spike as the value of Bi increases. However, note that the Biot number is specifically correlated to the coefficient of heat transfer *hf* ; therefore, it is conversely related to the thermal resistance of the current problem. Consequently, the heat resistance decreases as the Biot number increases, thereby increasing the heat transfer rate at the stretching/shrinking sheet and decreasing the temperature distribution (see Figure 15).
