Hybrid Nanoﬂuid Slip Flow over an Exponentially Stretching/Shrinking Permeable Sheet with Heat Generation

: An investigation has been done on the hybrid nanoﬂuid slip ﬂow in the existence of heat generation over an exponentially stretching/shrinking permeable sheet. Hybridization of alumina and copper with water as the base ﬂuid is considered. The mathematical model is simpliﬁed through the similarity transformation. A numerical solver named bvp4c in Matlab software is utilized to facilitate the problem-solving process and dual solutions are attained. The inﬂuences of several pertinent parameters on the main physical quantities of interest and the proﬁles are scrutinized and presented in the form of graphs. Through the stability analysis, only the ﬁrst solution is considered as the physical solution. As such, the ﬁndings conclude that the upsurges of volume fraction on the copper nanoparticle could enhance the skin friction coefﬁcient and the local Nusselt number.


Introduction
Heat transfer becomes the most important process with the evolution of advanced technology. Applications and industries involving heat generation, such as manufacturing, thermal power generation, transportation, chemical processes and many others, necessitate an efficient heat transfer performance to produce the best outcome. Nanofluid applications have been thoroughly reviewed by Sidik et al. [1,2] and also by Huminic and Huminic [3]. Historically, various investigations have been conducted to improvise the heat transfer capability of the transferring medium which is perceived as a heat transfer fluid.
The idea of dispersing high-thermal conductivity solid particles of microsize into the conventional fluid was first pioneered by Maxwell [4] and then continued by Hamilton and Crosser [5] to intensify the fluid thermal conductivity. Nevertheless, some limitations and flaws still occurred, such as the coagulation in the fluid flow passage. Owing to that limitation, a rapid investigation is being made and nanofluid is established. Choi and Eastman [6] invented this novel heat transfer fluid by dispersing the nanometer-sized solid particles into the base fluid, and it is believed these could overcome the coagulation of the flow passage due to its special feature.
However, along with the development of technologies, another class of heat transfer fluid named hybrid nanofluid is established as the extension to nanofluid that is invented by the dispersion of dual or multiple kinds of nanosized solid particles with a good thermal conductivity into a base fluid. Looking through the reviews of the previous investigations by the researchers and scientists [7][8][9][10][11], hybrid nanofluid is claimed to have a preferable capability in terms of thermal conductivity and heat transfer due to its synergistic effects, flow, which was then continued by Hayat et al. [48], who analyzed such effects on the rotating hybrid nanofluid. Eid and Nafe [42] also considered the presence of heat generation together with the slip effect in their research on hybrid nanofluid. Most of them indicated that the amplification of heat generation increases the temperature profile.
Therefore, in this present study, we intended to scrutinize the slip flow and heat transfer of a hybrid nanofluid of alumina and copper (hybrid nanoparticles) with water (base fluid) over an exponentially permeable stretching/shrinking sheet with the consideration of heat generation. Both velocity and thermal slips are being considered in the model. We believe the model of a hybrid nanofluid with such configurated effects and surface is still not being explored by the researchers, judging by the previous literature. The influences of various parameters over the main physical quantities are provided in the figure of the graph, and are discussed. An analysis of flow stability is also conducted, as the non-unique solutions are visible and the least eigenvalues that imply the stability of the flow are tabulated.

Mathematical Formulation
A steady boundary layer of hybrid nanofluid slip flow induced by an exponentially permeable stretching/shrinking sheet with heat generation is scrutinized. According to the physical model in Figure 1, u and v are the elements of velocity in x and y directions, respectively. The surface velocity is given as u w (x) = ce x/L and the wall mass transfer velocity is specified as v w (x) = v 0 e x/2L (v 0 < 0 is the mass suction and v 0 > 0 is the mass injection), λ > 0 is the stretching constant, λ < 0 is the shrinking constant,λ = 0 refers to a motionless sheet, T is temperature, T w = T ∞ + T 0 e x/2L is the sheet of varying temperature with constant T 0 and q = q 0 e x/L is the heat generation rate constant. heat generation. Murugesan and Kumar [46] also elucidated the same effect but with the flow of nanofluid over the exponentially stretching sheet. With the evolution of nanofluid to hybrid nanofluid, this heat generation/absorption effect has been considered by Hayat and Nadeem [47] with regard to the model of hybrid nanofluid with the three-dimensional flow, which was then continued by Hayat et al. [48], who analyzed such effects on the rotating hybrid nanofluid. Eid and Nafe [42] also considered the presence of heat generation together with the slip effect in their research on hybrid nanofluid. Most of them indicated that the amplification of heat generation increases the temperature profile.
Therefore, in this present study, we intended to scrutinize the slip flow and heat transfer of a hybrid nanofluid of alumina and copper (hybrid nanoparticles) with water (base fluid) over an exponentially permeable stretching/shrinking sheet with the consideration of heat generation. Both velocity and thermal slips are being considered in the model. We believe the model of a hybrid nanofluid with such configurated effects and surface is still not being explored by the researchers, judging by the previous literature. The influences of various parameters over the main physical quantities are provided in the figure of the graph, and are discussed. An analysis of flow stability is also conducted, as the non-unique solutions are visible and the least eigenvalues that imply the stability of the flow are tabulated.

Mathematical Formulation
A steady boundary layer of hybrid nanofluid slip flow induced by an exponentially permeable stretching/shrinking sheet with heat generation is scrutinized. According to the physical model in Figure 1, u and v are the elements of velocity in x and y directions, respectively. The surface velocity is given as  Therefore, concerning the mentioned boundary layer assumptions, the equations of continuity, momentum and energy are formulated as (Murugesan and Kumar [46], Ghosh and Mukhopadhyay [49], Waini et al. [16], Devi and Devi [14], Eid and Nafe [42]) Therefore, concerning the mentioned boundary layer assumptions, the equations of continuity, momentum and energy are formulated as (Murugesan and Kumar [46], Ghosh and Mukhopadhyay [49], Waini et al. [16], Devi and Devi [14], Eid and Nafe [42]) ∂u ∂x due to the boundary conditions (Ghosh and Mukhopadhyay [49], Mukhopadhyay and Andersson [50]) in which A 1 and B 1 are the slip factors for velocity and thermal, respectively. The correlation of the physical properties (Tayebi and Chamkha [51], Takabi and Salehi [52], Babu et al. [11], Ghalambaz et al. [53]) and the hybrid nanofluid's thermo-physical properties (Oztop and Abu-Nada [54]) is depicted in Tables 1 and 2, respectively. Here, the hnf, f, s1 and s2 subscripts refer to the hybrid nanofluid, base fluid and solid nanoparticle for the alumina and solid nanoparticle for copper, accordingly, and φ s1 and φ s2 symbolize the alumina and copper volume fraction parameters, separately. Table 1. Hybrid nanofluid physical properties.

Properties Hybrid Nanofluid Correlations
Density Heat Capacity Thermal Conductivity As to reduce the aforementioned governing equations, the similarity variables are employed (Waini et al. [16], Eid and Nafe [42]) which consequently gratify Equation (1) and transform Equations (2) and (3) into the subsequent equations Pr due to the transformed boundary conditions Mathematics 2021, 9, 30 5 of 20 such that Pr = µC p f /k f is the Prandtl number, β = 2q 0 L/c ρC p f is the heat generation parameter, A = A 1 µ hn f /ρ hn f e x/2L c/2ν f L is the dimensionless velocity slip parameter and B = B 1 e x/2L c/2ν f L is the dimensionless thermal slip parameter, and is the wall mass transfer parameter in which S > 0 (ν 0 < 0) equates to mass suction and S < 0 (ν 0 > 0) equates to mass injection.
The main physical quantities of interest (local skin friction coefficient C f and the local Nusselt number Nu x ) are mathematically computed as (Waini et al. [16], Eid and Nafe [42]) or by introducing the transform, we have where Re x = 2Lu w /ν f is the local Reynolds number.

Flow Stability
As it is attainable to compute dual solutions, the flow stability analysis is executed. The problem is regarded as unsteady in this analysis where the continuity equation (Equation (1)) is held, and the time derivative is introduced to Equations (2) and (3) u ∂u ∂x + v ∂u ∂y where the variable t here denotes time. Referencing the work reported by Merkin [55], Weidman et al. [56], Waini et al. [16] and Yan et al. [38], the dimensionless time-dependent variable transformation for Equations (11) and (12) is As such, with the utilization of the transformation, the subsequent equations are procured 1 Pr restricted to the conditions Mathematics 2021, 9, 30 6 of 20 The succeeding perturbation equation is introduced for the stability analysis of the similarity solutions, f (η) = f 0 (η) and θ(η) = θ 0 (η), where γ is an unidentified eigenvalue, while F(η) and G(η) are small relative to f 0 (η) and θ 0 (η), correspondingly (Weidman et al. [56], Waini et al. [16], Yan et al. [38]). Substituting Equation (17) into Equations (14)- (16) and setting the value of τ to zero, the initial decay or growth of disturbance in Equation (17) can be identified and the subsequent linearized equations are attained 1 Pr subject to the linearized conditions Since most of the linearized conditions (Equation (20)) are equal to zero, it is essential for F (∞) → 0 to relax and be replaced with F (0) = 1 to produce a realizable range of eigenvalues. Thus, to attain the unknown eigenvalues γ which act as the indicator to ascertain the stability of the flow, the bvp4c solver is utilized. The flow is contemplated to be stable when γ > 0, and vice versa.

Result and Discussions
The solutions for Equations (6) and (7) conditioned to Equation (8) are numerically enumerated by employing the bvp4c solver in Matlab software. The suitable preliminary guesses, preferable boundary layer thickness and various values of parameters need to be selected and well-adjusted in the coding function of the solver to compute the most precise results. Stretching (λ > 0) and shrinking (λ < 0) surfaces are considered in this study. The Prandtl number and nanoparticle volume fraction for alumina are specified to be fixed throughout this study as Pr = 6.2 and φ s1 = 0.01, respectively, while the other parameters, such as velocity A and thermal B slips, nanoparticle volume fraction for copper φ s2 , heat generation β and suction S, are set to be varied in order to study their impacts with regard to the boundary layer flow and heat transfer.
For verification purposes, the results presented in this study are contrasted with the past results revealed by Ghosh and Mukhopadhyay [49], Hafidzuddin et al. [57] and Waini et al. [16]. The numerical values of f (0) and −θ (0) are computed and compared in Table 3 for the case of viscous fluid (φ s1 = φ s2 = 0) and shrinking surface (λ = −1) with S = 3 and Pr = 0.7. We observe that the present and previous results are found to be in a reasonable correlation. Thus, we guarantee that the technique, as well as the results provided in this present study, are valid and acceptable. Table 3. Tabulation of f (0) and −θ (0) for viscous fluid (φ s1 = φ s2 = 0) when Pr = 0.7, S = 3 and A = B = β = 0 for the shrinking surface (λ = −1) case. The non-unique or the dual solutions of Equations (6) and (7) are perceived to exist in this study for certain values of parameters, as can be observed in Figures 2-9 for the plots of skin friction coefficient Re 1/2 x C f and the local Nusselt number Re −1/2 x Nu x , and also in                               In Figure 2, the skin friction coefficient Re 1/2 x C f is noted to be in a decreasing trend as the velocity slip parameter A escalates for the first solution in a certain range (−1.4419 ≤ λ < 0) in the shrinking surface (λ < 0) region, but the opposite trend occurs when approaching the stretching surface (λ > 0) region and onwards. However, for the local Nusselt number Re −1/2 x Nu x , it is noticed from Figure 3 that the enhancement in the value of the velocity slip parameter A promotes the increment in the local Nusselt number Re −1/2 x Nu x for the first solution, but when approaching the stretching surface (λ > 0) region, the opposite trend occurs. Meanwhile, for the second solution of the skin friction coefficient Re 1/2 x C f and the local Nusselt number Re −1/2 x Nu x , a similar trend is also noticed to occur like the first solution, but a change of trend happens within a certain range of shrinking surface (λ < 0) region and the stretching surface region (λ > 0) onwards, unlike in the first solution. The skin friction coefficient Re 1/2 x C f is also perceived to increase as the value of the stretching/shrinking parameter λ diminishes, while the obverse occurs for the local Nusselt number Re −1/2 x Nu x for both the first and second solutions. x Nu x to decrease, but it is seen to increase as the stretching/shrinking parameter rises λ for both the first and second solutions.
Moreover, for the influence of copper volume fraction φ s2 (= 0.001, 0.005, 0.01) over the skin friction coefficient Re 1/2 x C f , and the local Nusselt number Re −1/2 x Nu x against the stretching/shrinking parameter λ when S = 2.4, A = B = 0.2, β = 0.02, is shown in Figures 6 and 7, and the same against the suction parameter S when λ = −1, A = B = 0.2, β = 0.02 is shown in Figures 8 and 9. It is noted that in Figures 6 and 7, the existence of dual solutions is visible for φ s2 (= 0.001, 0.005, 0.01) when λ c1 = −1.3686, λ c2 = −1.4016, λ c3 = −1.4419, which indicates that the bifurcation point of the boundary layer eventuates at the shrinking surface (λ < 0) region and the presence of the copper volume fraction delays the bifurcation of the boundary layer.
The upsurge in the copper volume fraction parameter φ s2 leads to the increase in the skin friction coefficient Re 1/2 x C f at the shrinking surface (λ < 0) region for the first solution, but leads to a reduction for the second solution. We also notice a change in the trend of the skin friction coefficient Re 1/2 x C f for the first solution towards the copper volume fraction parameter φ s2 when λ > 0. Meanwhile, the plot of the local Nusselt number Re −1/2 x Nu x is observed to rise when there is an upsurge in the copper volume fraction parameter φ s2 for the first and second solutions. The plots of the skin friction coefficient Re 1/2 x C f and local Nusselt number Re −1/2 x Nu x are noticed to increase and decrease, respectively, when the stretching/shrinking parameter λ decreases for both solutions.
In Figures 8 and 9, dual solutions are perceived to exist for φ s2 (= 0.001, 0.005, 0.01) when S c1 = 2.0775, S c2 = 2.0551, S c3 = 2.0290, which indicates that only certain appropriate values of suction are required to achieve dual solutions for the shrinking surface. In these figures, the increases in copper volume fraction φ s2 and suction S parameters augment the skin friction coefficient Re 1/2 x C f for the first solution but reduce the skin friction coefficient Re 1/2 x C f for the second solution. Concurrently, the local Nusselt number Re −1/2 x Nu x is noted to augment when the copper volume fraction φ s2 and suction S parameters increase for both of the solutions. From all of the figures (Figures 2-9), we noted that the plots of the skin friction coefficient and the local Nusselt number against both the stretching/shrinking and the suction parameters are all moving towards the left side of the domain, either up or down on the left side, to achieve the critical point.
The velocity f (η) and temperature θ(η) profiles are also provided as in Figures 10-17. The increment of velocity slip A (= 0.2, 0.4, 0.6) and suction S (= 2.2, 2.3, 2.4) parameters causes the velocity profile f (η) to increase for the first solution, but it reduces for the second solution and for both solutions of the temperature profile θ(η) as depicted in Figure 10, Figure 11, Figure 12, Figure 13, separately. For the first realizable solution, physically, the presence of velocity slip at the boundary of the shrinking surface reduces the flow resistance which causes more flow to slip through the surface and increase the flow velocity, thus a high degree of velocity slip is required to speed up the flow velocity and vice versa. For the copper volume fraction parameter φ s2 (= 0.001, 0.005, 0.01), the upsurge in value amplifies both the velocity f (η) and temperature θ(η) profiles for both of the solutions, as exemplified in Figures 14 and 15, but in the velocity profile f (η) for the second solution, dissimilar behavior is noticed at a certain range. In Figures 16 and 17, the influences of the thermal slip parameter B (= 0.2, 0.4, 0.6) and heat generation β (= 0.02, 0.04, 0.06) parameter over the temperature profile θ(η) are portrayed. Since these parameters do not affect the velocity profile f (η), it is only necessary to provide the temperature profile θ(η) in this present study. The temperature profile θ(η) is noticed to increase when the thermal slip parameter B decreases and when the heat generation parameter β increases. The appropriate amounts of thermal slip and heat generation parameters need to be controlled to obtain the required heat transfer performance, as both parameters produce a different impact that could offset each other.
As dual solutions exist, it is essential to execute the stability analysis of the flow. The stability analysis is conducted as has been discussed in the previous section. The computations of the lowest eigenvalues γ for A = 0.2 and A = 0.6 for different λ when S = 2.4, φ s2 = 0.01, B = 0.2 and β = 0.02 are tabulated in Table 4. It is realized that the first solution gives the positive value of eigenvalues γ 1 , meanwhile the second solution gives the negative value of eigenvalues γ 2 . As according to Merkin [55] and Weidman et al. [56], the positive eigenvalues γ 1 imply that the flow is real and stable due to the initial decay of disturbance as the time passes, and the negative eigenvalues γ 2 imply that the flow is unreal and unstable due to the initial growth of disturbance. It is also noticeable from the table that as the value of λ → λ c , the value of γ → 0 . So, as supported by the statement above, it can be affirmed that the flow is stable for the first solution whilst it is unstable for the second solution.

Conclusions
The slip flow of a hybrid nanofluid over a permeable exponentially stretching/shrinking sheet with the existence of heat generation has been numerically scrutinized by using the bvp4c solver. The present findings can be summarized as follows: 1.
Dual solutions are achievable at the shrinking surface region as well as with the appropriate amount of suction, and only the first solution is stable.

2.
The increment of the velocity slip parameter depreciates the skin friction coefficient but enhances the local Nusselt number in the shrinking surface region (first solution).

3.
The skin friction coefficient is intensified (first solution) when the values of the copper volume fraction and suction parameters increase in the shrinking surface region.

4.
The local Nusselt number is augmented (first and second solutions) when the copper volume fraction and suction parameters increase, but is diminished (first and second solutions) as the thermal slip and heat generation parameters increase.

5.
The rise in the stretching/shrinking parameter reduces the skin friction coefficient but augments the local Nusselt number. 6.
The augmentation of the velocity slip, suction and copper volume fraction parameters augments the velocity profile (first solution). 7.
The temperature profile can be enhanced (first and second solutions) with the enhancement of the copper volume fraction and heat generation parameters, but can be diminished (first and second solutions) with the increment of the velocity and thermal slips and suction parameters.