Dual Solutions and Stability Analysis of Magnetized Hybrid Nanoﬂuid with Joule Heating and Multiple Slip Conditions

: This paper investigates the steady, two dimensional, and magnetohydrodynamic ﬂow of copper and alumina / water hybrid nanoﬂuid on a permeable exponentially shrinking surface in the presence of Joule heating, velocity slip, and thermal slip parameters. Adopting the model of Tiwari and Das, the mathematical formulation of governing partial di ﬀ erential equations was constructed, which was then transformed into the equivalent system of non-linear ordinary di ﬀ erential equations by employing exponential similarity transformation variables. The resultant system was solved numerically using the BVP4C solver in the MATLAB software. For validation purposes, the obtained numerical results were compared graphically with those in previous studies, and found to be in good agreement, as the critical points are the same up to three decimal points. Based on the numerical results, it was revealed that dual solutions exist within speciﬁc ranges of the suction and magnetic parameters. Stability analysis was performed on both solutions in order to determine which solution(s) is / are stable. The analysis indicated that only the ﬁrst solution is stable. Furthermore, it was also found that the temperature increases in both solutions when the magnetic parameter and Eckert number are increased, while it reduces as the thermal slip parameter rises. Furthermore, the coe ﬃ cient of skin friction and the heat transfer rate increase for the ﬁrst solution when the magnetic and the suction parameters are increased. Meanwhile, no change is noticed in the boundary layer separation for the various values of the Eckert number in the heat transfer rate.


Introduction
Previous studies have proven that fluids play a vital role in enhancing the heat transfer rate in many engineering and industrial systems. It is also noticed from the previous literature that regular fluids, such as water, ethylene glycol, etc., must keep thermal conductivity low in order to improve the heat transfer rate. Many researchers have attempted to find new fluids that can further enhance the heat transfer rate. Choi and Eastman [1] were the first that introduced another kind of fluid known as "nanofluid" and claimed that the heat transfer rate of nanofluid outperformed that of regular fluids. Nanofluid can be defined as a mixture of nano-sized particles of solid in a base fluid. Later, many researchers worked on nanofluids experimentally and theoretically [2][3][4][5][6][7][8]. Consequently, these investigations lead to greater efforts to enhance the heat transfer rate, yet nobody can conclude and claim that a mixture of a particular type of base fluid with a particular kind of nanoparticle can produce the highest rate of heat transfer [9][10][11][12][13][14]. Continuous attempts have been made to improve the rate of heat transfer using various strategies. One of the approaches is mixing two different types of solid nanoparticle in a solitary base fluid, which produces a special kind of nanofluid known as hybrid nanofluid. This idea was initiated by Suresh et al. [15,16] by considering alumina and copper as the solid particles with a base fluid of water. Furthermore, they claimed that the rate of heat transfer of hybrid nanofluid is higher at the surface compared to that of simple fluid and nanofluid. Their discovery is a new research direction in fluid dynamics. Since then, numerous researchers have been focusing on hybrid nanofluid in their work and have discovered some interesting and promising results, as can be seen in these articles [17][18][19][20][21][22]. In 2016, Devi and Devi [23] contemplated the copper-alumina/water hybrid nanofluid numerically. They compared the numerical results of thermal conductivity with the experimental work of Suresh et al. [16] and found an astounding agreement between them. In the same year, the 3D magnetohydrodynamics MHD flow of hybrid nanofluid was considered by Devi and Devi [24], and their numerical results coincided with the experimental results obtained by Suresh et al. [15]. Because of these encouraging findings, the newly found thermophysical properties have been considered by many researchers in their studies, such as Das et al. [25], Mehryan et al. [26], Chamkha et al. [27], Hayat et al. [28], Saba et al. [29], Kumar and Sarkar [30], Kassai [31], and Ghalambaz et al. [32]. Recently, Waini et al. [33] considered the unsteady flow of hybrid nanofluid, by mixing copper and alumina (nanoparticles) with water (base fluid). They noticed that double solutions appeared in specific ranges of the unsteadiness parameter. Lund et al. [34] examined the steady flow of a hybrid nanofluid in the existence of a dissipation function effect and found double solutions. In addition, stability analysis on the solutions was done and demonstrated that only the first solution is stable. Meanwhile, the axisymmetric flow of electrically conducting hybrid nanofluid was inspected by Khashi'ie et al. [35], in whose work the occurrence of dual solutions is noticed within certain ranges of physical parameters. The stability analysis also specified that only the first solution is stable. Anuar et al. [36] showed the existence of dual solutions in the stagnation point-flow of water-based hybrid nanofluid with the velocity slip effect and concluded that a non-uniqueness of solutions exists over the exponentially shrinking surface. Their work was then continued by Waini et al. [37] by considering the effects of MHD and thermal radiation without stagnation point flow and slip effects. It was noticed that double solutions occur over the shrinking and stretching surfaces within specific ranges of the mass suction parameter.
Nowadays, fluid flow over the shrinking surface has gotten the attention of researchers who are intrigued to find the multiple solutions. There are greater possibilities for the existence of dual solutions when fluid is flowing on shrinking surfaces [38]. It appears that Miklavčič and Wang [39] are the main researchers that considered the flow of viscous fluid over the shrinking surface and reasoned that flow over the shrinking surface is not going to exist due to the unconfined vorticity. In order to keep the flow, sufficient suction on the surface is required. Afterwards, Fang et al. [40] examined viscous fluid flow over the shrinking sheet including the second slip parameter, and revealed that the impact of drag force can be decreased by applying higher mass suction on the surface. Meanwhile, Jahan et al. [41] analyzed the unsteady flow of nanofluid over the shrinking sheet and discovered double solutions.
Hybrid nanofluid flow on the stretching/shrinking surface, with the effects of transpiration, was explored by Waini et al. [42]. It is also noticed from the previous literature that regular fluids, such as water, ethylene glycol, etc., must keep thermal conductivity low in order to improve the heat transfer rate [43,44]. Some current developments on nanofluids and hybrid nanofluids over the shrinking surface for multiple solutions can be seen in these articles [43][44][45][46][47].
A huge overview of the published work demonstrates that the effect of slip conditions on nanofluid flow has not been given a lot of consideration, especially on hybrid nanofluids. Numerous significant applications of fluids show boundary slip conditions, such as the perfecting of heart valves and inside cavities, and the cleaning of artificial heart valves [48]. It is worthwhile to state that the no-slip condition is not always valid in reality. In simple words, the velocity of the slip condition can be explained, as the fluid would not have zero velocity with regards to contact with the solid boundary. Similarly, a thermal slip condition can be explained. Andersson [49] is probably the first who introduced the concept of the slip effect on the boundary layer flow. Uddin et al. [50] examined nanofluid with the impacts of slip conditions and Darcian porous medium, and inferred that the temperature slip parameter has an inverse relationship with the temperature and velocity profiles. The MHD flow of a hybrid nanofluid in the presence of the slip condition was considered by Nadeem and Abbas [51], and they adopted the thermophysical properties of Devi and Devi [23,24] to solve the system of ordinary differential equations (ODEs). Using the same properties, Iftikhar et al. [52] also investigated the steady MHD flow of hybrid nanofluid with the slip condition. Based on published literature, it can be concluded that the thermophysical properties of Devi and Devi have been used broadly, as these properties match with the experimental results. Taking advantage of this situation, in this study, we also considered the thermophysical properties of Devi and Devi [23,24] in dealing with hybrid nanofluid, and anticipated that our results would help in understanding hybrid nanofluid effectively without having to conduct costly experimental studies. This paper has three main objectives. First, to extend the work of Anuar et al. [36] and Waini et al. [37] by incorporating Joule heating effects with velocity and thermal slip conditions. By including the Joule heating effect, it's very hard to find the solution due to the non-linearity in the temperature equations. Therefore, many researchers did not consider it in their study. Second, to obtain the maximum number of multiple solutions due to the existence of non-linearity. Third, to find a stable solution by performing stability analysis. According to the best of our knowledge, no such study has been contemplated previously. A 10% of Al 2 O 3 volume fraction with a 0.1% ≤ Cu ≤ 10% volume fraction was chosen for the numerical computation, and produced interesting results. It is expected that the current studies will provide good benefits to the researchers who are working on the hybrid nanofluid experimentally, and it is also expected that these results will reduce the cost of the experimental work in the future.

Mathematical Formulation
There have been considered the incompressible 2D, MHD, and steady boundary layer flow of hybrid nanofluid with the effects of Joule heating, velocity, and thermal slip conditions over the exponentially shrinking surface (see Figure 1) without the viscous dissipation effect. By considering all assumptions, the governing mass, momentum, and energy conservations can be expressed as [35,37]: Here, u and v are the corresponding velocities of the x-axis and y-axis; T is the temperature of fluid; T w (x) is the temperature of the surface and is explained as T w (x) = T ∞ + T 0 e    The following variables of similarity transformation are used to reduce the system into ODEs [37]

Properties Hybrid Nanofluid
Dynamic viscosity Density Thermal conductivity Electrical conductivity Along with the boundary conditions where primes stand for the derivatives with respect to η, 2l is the velocity slip parameter, and δ T = D 1 U w 2ϑ f l is the thermal slip parameter.
The important physical quantities of interest are the skin friction coefficient, C f , and local Nusselt number, Nu x , explained as Processes 2020, 8, 332 6 of 17 Applying similarity variables (5) yields

Stability Analysis
In order to perform the stability analysis of solutions of Equations (6) and (7) along with boundary conditions (9), the unsteady forms of the Equations (2) and (3) are required to test the features of the temporal stability analysis by introducing the new time-dependent dimensionless variable, τ = U w 2l e x l .t. In this regard, the steps of Merkin [55], Lund et al. [56], and Weidman et al. [57] are adopted. Thus, we get ∂u ∂t Equation (5) can take the following form with the new transformation variable of τ By using Equation (14) into Equations (12) and (13), we get Subject to boundary conditions To check the stability of steady flow solutions where f (η) = f 0 (η) and θ(η) = θ 0 (η) of satisfying the boundary value problem (6-9), we have where ε is the unknown eigenvalue; during the process of finding of ε, it will generate the unlimited set of the eigenvalues, ε 1 < ε 2 < ε 3 . . . . Furthermore, F(η, τ) and G(η, τ) are the small related values of θ 0 (η) and f 0 (η), where f (η) = f 0 (η) and θ(η) = θ 0 (η) can be obtained by setting τ = 0. Basically, these show the solutions of Equations (6) and (7). Accordingly, the flowing linear eigenvalue problems need to be solved. Subject to the reduced boundary conditions To get the smallest eigenvalue, ε 1 , one of them, F 0 (η) → 0, G 0 (η) → 0 as η → ∞ , is required to relax into an initial boundary condition, as suggested by Haris et al. [58]. In this problem, F 0 (η) → 0 as η → ∞ is relaxed to F 0 (0) = 1.

Results and Discussion
The system of non-linear ODEs (6-7) along boundary conditions (9) was solved by utilizing BVP4C in MATLAB programming. In order to get the solutions with good precision, we had to make some initial guesses at the initial mesh points by the changing of step-sizes. In the whole study, the tolerance was set at 10 −6 in order to get good accuracy in the solutions. Waini et al. [33] and Lund et al. [54] described this method in detail. As there exist two solutions, accordingly, more initial guesses were required to acquire the solution subject to fulfill boundary conditions asymptotically at η → ∞ . In the entire study, η → ∞ was selected as η = 8 for both solutions.
The numerical calculations were carried out for the various physical included parameters. In Figure 2, a graphical comparison of √ ReC f with Waini et al. [37] has been made for the validation of our numerical method, and shows an excellent agreement. The behavior of the graph and the critical values of the suction parameter, S c , are equivalent up to three decimal places, as referenced in the paper of Waini et al. [37] (alluding to Figure 2 of [37]). Henceforth, the present method can certainly be used, and the generated results are reliable and correct. Additionally, the numerical values of f (0) and −θ (0) for the different values of the magnetic parameter (M), the solid volume fraction of copper parameter (φ Cu ), the Prandtl (Pr) and Eckert (Ec) numbers, the Suction parameter (S), the velocity slip parameter (δ), and the temperature slip parameter (δ T ) are given in Tables 3 and 4.

Results and Discussion
The system of non-linear ODEs (6-7) along boundary conditions (9) was solved by utilizing BVP4C in MATLAB programming. In order to get the solutions with good precision, we had to make some initial guesses at the initial mesh points by the changing of step-sizes. In the whole study, the tolerance was set at 10 in order to get good accuracy in the solutions. Waini et al. [33] and Lund et al. [54] described this method in detail. As there exist two solutions, accordingly, more initial guesses were required to acquire the solution subject to fulfill boundary conditions asymptotically at → ∞.
In the entire study, → ∞ was selected as  = 8 for both solutions.
The numerical calculations were carried out for the various physical included parameters. In Figure 2, a graphical comparison of √ with Waini et al. [37] has been made for the validation of our numerical method, and shows an excellent agreement. The behavior of the graph and the critical values of the suction parameter, , are equivalent up to three decimal places, as referenced in the paper of Waini et al. [37] (alluding to Figure 2 of [37]). Henceforth, the present method can certainly be used, and the generated results are reliable and correct. Additionally    Physically, it is due to the fact that the higher impact of the magnetic number creates resistance in the fluid flow, and vorticity is smothered, as seen in the first solution. The skin friction coefficient, f (0), and rate of heat transfer, −θ (0), increase when the magnetic number, M, increases for the first solution. These increments are due to the fact that the magnetic field creates the Lorentz force. On the other hand, they decline in the second solution.  Figures 5 and 6 show the variation of skin friction coefficient, f (0), and heat transfer rate, −θ (0), against the suction parameter, S, for fixed values of the solid volume fraction of copper, φ Cu . When φ Cu increases, the critical values of suction, S c , increment towards the left, which supports the delay in the boundary layer separation. The critical values of S for φ Cu = 0, 0.04 and 0.1 are S c1 = 1.7373, S c2 = 1.6796, and S c3 = 1.6312 respectively. Furthermore, the presence of dual solutions is conceivable when S ≥ S ci where i = 1, 2, 3, whereas no solution exits when S < S ci . It is observed that for the fixed values of S, the skin friction coefficient, f (0), increases with the expanding of φ Cu for the first solution, whereas the inverse trend is seen in the second solution. Moreover, the heat transfer rate, −θ (0), increases when the suction parameter is increased, by keeping the fixed values of φ Cu in the first solution. It can be explained as so: suction produces the drag force, which slows down the fluid flow, and as a result, the heat transfer rate increases. Meanwhile, the rate of heat transfer diminishes gradually in the second solution.
parameter ( ) are illustrated in Figures 3 and 4 separately. It is observed from these figures that dual solutions exist for Equations (6) and (7) along the boundary condition (9) in the specific ranges of . The respective critical values for = 0.001, 0.04 and 0.1 are = 0.4847, = 0.4112, and = 0.3261 where solutions exist. There exist two ranges of solution which rely on the values of , no solution exists when < where = 1, 2, 3 , and dual solutions exist when . Furthermore, it is shown that decays with the increasing values of , which show the delay in the boundary layer separation. Physically, it is due to the fact that the higher impact of the magnetic number creates resistance in the fluid flow, and vorticity is smothered, as seen in the first solution. The skin friction coefficient, ′′(0), and rate of heat transfer, − ′(0), increase when the magnetic number, , increases for the first solution. These increments are due to the fact that the magnetic field creates the Lorentz force. On the other hand, they decline in the second solution.    The effect of the velocity slip parameter, δ, on the skin friction coefficient, f (0), for numerous values of S is demonstrated in Figure 7. It is observed that the skin friction coefficient, f (0), increases initially, and then decreases for the higher values of δ in the first solution; practically, it shows that the viscosity of the hybrid nanofluid enhances initially, and then starts to decrease after the intensive impact of the S parameter in the boundary layer. However, a contradictory nature of the skin friction coefficient, f (0), is noticed for the second solution. Moreover, the critical values of S for the higher values of δ are S c1 = 1.5815, S c2 = 1.6686, and S c3 = 1.7648 where solutions exist, while no solution exists beyond these critical values. The effect of the Eckert number, Ec, on the heat transfer rate, −θ (0), is portrayed in Figure 8. The reduction of heat transfer is noticed in the two solutions as Ec is enhanced. However, there is no effect of the higher values of Ec on the delaying boundary layer separation, as the critical values of S are the same. Dual solutions exist when S ≥ S c = 1.6686 and no solution exists when S < S c . It is worthwhile to mention that the effect of Ec on f (0) is not important, as we have gotten the same values of f (0) for the higher values of Ec. Figure 9 displays the variation of the heat transfer rate, −θ (0), for numerous values of δ T . The same behavior of −θ (0) is observed as it is demonstrated in Figure 8.  Figures 5 and 6 show the variation of skin friction coefficient, ′′(0), and heat transfer rate, − ′(0), against the suction parameter, , for fixed values of the solid volume fraction of copper, . When increases, the critical values of suction, , increment towards the left, which supports the delay in the boundary layer separation. The critical values of for = 0, 0.04 and 0.1 are = 1.7373, = 1.6796, and = 1.6312 respectively. Furthermore, the presence of dual solutions is conceivable when where = 1, 2, 3, whereas no solution exits when < . It is observed that for the fixed values of , the skin friction coefficient, ′′(0), increases with the expanding of for the first solution, whereas the inverse trend is seen in the second solution. Moreover, the heat transfer rate, − ′(0), increases when the suction parameter is increased, by keeping the fixed values of in the first solution. It can be explained as so: suction produces the drag force, which slows down the fluid flow, and as a result, the heat transfer rate increases. Meanwhile, the rate of heat transfer diminishes gradually in the second solution.   Figure 8. The reduction of heat transfer is noticed in the two solutions as is enhanced. However, there is no effect of the higher values of on the delaying boundary layer enhanced. However, there is no effect of the higher values of on the delaying boundary layer separation, as the critical values of are the same. Dual solutions exist when = 1.6686 and no solution exists when < . It is worthwhile to mention that the effect of on ′′(0) is not important, as we have gotten the same values of ′′(0) for the higher values of . Figure 9 displays the variation of the heat transfer rate, − ′(0), for numerous values of . The same behavior of − ′(0) is observed as it is demonstrated in Figure 8.    Figures 10 and 11 show the effect of φ Cu and φ Al 2 O 3 on the velocity profile, f (η), and temperature profile, θ(η). It is found that the velocity and thickness of the momentum boundary layer decrease for both solutions with increments in φ Cu when φ Al 2 O 3 = 0.1 is kept as a constant. In contrast, the momentum boundary layer becomes thicker in the subsequent solution when φ Al 2 O 3 is increased for the fixed value of φ Cu = 0.1, while a reverse trend of velocity profile is noticed for the first solution. Furthermore, the thickness of the thermal boundary layer and the temperature distributions increase in both solutions with rising values of φ Cu and φ Al 2 O 3 . It is examined that by varying φ Cu , values span a greater range of temperatures in the fluid flow from the surface as compared to φ Al 2 O 3 ; this behaviour of the temperature profiles supports our assumptions of keeping φ Al 2 O 3 fixed in the whole study. The thermal boundary layer thickness develops as the Eckert and magnetic numbers are increased in both solutions, as seen in Figure 12. This rising of the temperature is due to the higher kinetic energy, which is directly proportional to the Eckert number. Furthermore, the Lorentz force also causes the temperature of the fluid to spread from the surface towards the flow.  ; this behaviour of the temperature profiles supports our assumptions of keeping fixed in the whole study. The thermal boundary layer thickness develops as the Eckert and magnetic numbers are increased in both solutions, as seen in Figure 12. This rising of the temperature is due to the higher kinetic energy, which is directly proportional to the Eckert number. Furthermore, the Lorentz force also causes the temperature of the fluid to spread from the surface towards the flow.  In order to indicate the physically realizable solution, the implementation of the stability analysis is necessary to conduct in the study when a non-uniqueness of solutions exists. Normally, the first solution is referred to as the physical solution, as it satisfies the far-field boundary condition, but it cannot be said which solution is the physical solution without performing the stability analysis of solutions. The indicated solution might be the second solution. Therefore, the performing of stability analysis of solutions should be considered to stop wrong predictions regarding the flow characteristics and solutions. In this study, the system of the eigenvalue problem (19)(20) was solved with the help of BVP4C code in the MATLAB software in order to obtain the values of the smallest eigenvalue, ε 1 . The signs of the smallest eigenvalue, ε 1 , help to indicate a stable solution. If the sign of ε 1 is positive, it is implied that the flow is stable and showing an initial growth of decay, while if it is negative, it means that the flow is unstable and indicates the initial growth of disturbance. Positive values of ε 1 can be seen in the upper part of Figure 13, and negative values of ε 1 in the lower part. Therefore, the first solution is stable and the second is unstable. Furthermore, it is also noticed that ε 1 0 for both solutions when M = M c , which validates the current formulation of the problem and also proves that the real solution is the first solution.   In order to indicate the physically realizable solution, the implementation of the stability analysis is necessary to conduct in the study when a non-uniqueness of solutions exists. Normally, the first solution is referred to as the physical solution, as it satisfies the far-field boundary condition, but it cannot be said which solution is the physical solution without performing the stability analysis of solutions. The indicated solution might be the second solution. Therefore, the performing of stability analysis of solutions should be considered to stop wrong predictions regarding the flow characteristics and solutions. In this study, the system of the eigenvalue problem (19-20) was solved of is positive, it is implied that the flow is stable and showing an initial growth of decay, while if it is negative, it means that the flow is unstable and indicates the initial growth of disturbance. Positive values of can be seen in the upper part of Figure 13, and negative values of in the lower part. Therefore, the first solution is stable and the second is unstable. Furthermore, it is also noticed that ≅ 0 for both solutions when = , which validates the current formulation of the problem and also proves that the real solution is the first solution.

Conclusions
The steady MHD flow of hybrid nanofluid with the effects of Joule heating and slip conditions has been investigated over the shrinking surface. After the transformation of the similarity variables, equations were solved by employing BVP4C. The pointwise conclusions of the present study are as follows: 1.
There exist two ranges of solution, namely dual solutions and no solution.

2.
Dual solutions do not exist beyond the critical values (S c , M c ) of the parameters. 3.
The existence of dual solutions is possible in certain dimensions of the suction parameter S.

4.
Due to the effect of Joule heating, the dual solutions also depend on certain ranges of the magnetic parameter, M.

5.
The skin friction coefficient, f (0), enhances for the first solution when the S and M parameters are increased, while f (0) reduces for the higher effect of the velocity slip factor, δ.

6.
The heat transfer rate, −θ (0), reduces with increments in the Eckert number, Ec, and the thermal slip parameter, δ T ; however, Ec and δ T do not affect the boundary layer separation. 7.
The temperature and thermal boundary layer thickness have direct relationships with Ec for both solutions. 8.
Positive smallest eigenvalues indicate the initial decay of the disturbance, and that the flow becomes the stable flow. 9.
The stability analysis indicates that the real solution is the first solution.