Flow towards a Stagnation Region of a Vertical Plate in a Hybrid Nanoﬂuid: Assisting and Opposing Flows

: This study investigates a hybrid nanoﬂuid ﬂow towards a stagnation region of a vertical plate with radiation effects. The hybrid nanoﬂuid consists of copper (Cu) and alumina (Al 2 O 3 ) nanoparticles which are added into water to form Cu-Al 2 O 3 /water nanoﬂuid. The stagnation point ﬂow describes the ﬂuid motion in the stagnation region of a solid surface. In this study, both buoyancy assisting and opposing ﬂows are considered. The similarity equations are obtained using a similarity transformation and numerical results are obtained via the boundary value problem solver (bvp4c) in MATLAB software. Findings discovered that dual solutions exist for both opposing and assisting ﬂows. The heat transfer rate is intensiﬁed with the thermal radiation (49.63%) and the hybrid nanoparticles (32.37%).


Introduction
The phenomenon of the flow on a stagnation region commonly occurs in aerodynamic industries and engineering applications. To name a few, such applications are polymer extrusion, drawing of plastic sheets, and wire drawing. Hiemenz [1] was the first researcher to consider the boundary layer flow toward a stagnation point on a rigid surface. Besides this, the axisymmetric flow was considered by Homann [2], whereas the oblique stagnationpoint flow was studied by Chiam [3]. Further, Merkin [4] studied a similar problem by considering the mixed convection flow. He discovered that the solution is not unique for the opposing flow case. However, Ishak et al. [5] exposed that the dual solutions exist for both opposing and assisting flows, and these behaviours were also reported by several researchers [6][7][8][9].
In 1995, Choi and Eastman [10] presented a new type of heat transfer fluid called nanofluid, which is a mixture of single type nanoparticles and the base fluid, to enhance the thermal conductivity. Some works on such fluids can be found in [11][12][13][14][15][16]. Recently, some studies have shown that advanced nanofluids composed of other types of nanoparticles mixed with regular nanofluids could improve their thermal properties, and this mixture is termed "hybrid nanofluid". The earlier experimental works on the hybrid nanofluid have been done by Turcu et al. [17], Jana et al. [18], and Suresh et al. [19]. Besides, the numerical studies on the hybrid nanofluid flow were studied by Devi and Devi [20]. They observed that the heat transfer rate of the hybrid nanofluid is higher than that of the regular nanofluid. Moreover, the non-uniqueness of the solutions in the hybrid nanofluid flow was examined by Waini et al. [21][22][23][24][25][26][27] Other physical aspects were considered by several authors [28][29][30][31][32][33][34][35]. Furthermore, the review papers can be found in [36][37][38][39][40][41]. Different from the above-mentioned studies, this paper considers the assisting and opposing buoyant flows of a hybrid nanofluid containing Al 2 O 3 -Cu hybrid nanoparticles when the effect of thermal radiation is taken into consideration. The governing equations along with the boundary conditions are transformed into a system of ordinary differential equations using a similarity transformation. The system of equations is then solved numerically using the boundary value problem solver (bvp4c) in MATLAB software. Most importantly, in this study, two solutions are discovered for both opposing and assisting flows. Then, further analysis is performed to study the temporal stability of these solutions as time evolves.

Mathematical Formulation
Consider the flow configuration as shown in Figure 1. The free stream velocity is U(x) = ax and the surface temperature is T w (x) = T ∞ + bx, where a and b are constants. Meanwhile, the ambient temperature T ∞ is assumed to be constant. Accordingly, the hybrid nanofluid equations are as follows ( [5,14]): where u and v represent the velocity components along the xand yaxes. Besides, g and q r are the acceleration caused by the gravity and the radiative heat flux, respectively. Meanwhile, the temperature of the hybrid nanofluid is given by T.
Mathematics 2021, 9, x FOR PEER REVIEW 2 of 16 Different from the above-mentioned studies, this paper considers the assisting and opposing buoyant flows of a hybrid nanofluid containing Al2O3-Cu hybrid nanoparticles when the effect of thermal radiation is taken into consideration. The governing equations along with the boundary conditions are transformed into a system of ordinary differential equations using a similarity transformation. The system of equations is then solved numerically using the boundary value problem solver (bvp4c) in MATLAB software. Most importantly, in this study, two solutions are discovered for both opposing and assisting flows. Then, further analysis is performed to study the temporal stability of these solutions as time evolves.

Mathematical Formulation
Consider the flow configuration as shown in Figure 1. The free stream velocity is ( ) = and the surface temperature is ( ) = + , where and are constants. Meanwhile, the ambient temperature is assumed to be constant. Accordingly, the hybrid nanofluid equations are as follows ( [5,14]): subject to where and represent the velocity components along the -and -axes. Besides, and are the acceleration caused by the gravity and the radiative heat flux, respectively. Meanwhile, the temperature of the hybrid nanofluid is given by .  The expression of the radiative heat flux is ( [42,43]): where σ * and k * denote the Stefan-Boltzmann constant and the mean absorption coefficient, respectively. Following Rosseland [42], after employing a Taylor series, one gets Then, the Equation (3) turns to [43]: Further, the thermophysical properties can be referred to in Tables 1 and 2. Data from these tables are adapted from Oztop and Abu-Nada [13], Devi and Devi [20], and Waini et al. [21]. Note that ϕ 1 (Al 2 O 3 ) and ϕ 2 (Cu) are the nanoparticles volume fractions, and the subscripts n1 and n2 are corresponded to their solid components, while the subscripts f , n f , and hn f signify the base fluid, nanofluid, and hybrid nanofluid, respectively.
To get a similarity solution, we employ the following similarity transformation ( [5,14]): where ψ is the stream function defined as u = ∂ψ/∂y and v = − ∂ψ/∂x, then one gets  Table 2. Thermophysical properties of nanofluid and hybrid nanofluid.

Properties Nanofluid Hybrid Nanofluid
Dynamic viscosity Density Heat capacity Thermal conductivity Furthermore, the continuity equation, i.e., Equation (1), is identically satisfied. Now, Equations (2) and (6) respectively reduce to: subject to the boundary conditions: where ( ) represents the differentiation with respect to η, Pr is the Prandtl number, R and λ signify the radiation and the mixed convection parameters, given by:

corresponds to the local Grashof number and
Re x = ax 2 /ν f stands for the local Reynold's number. Note that λ < 0 signifies the opposing and λ > 0 signifies the assisting flows, while the forced convection flow (no buoyancy effects) is given by λ = 0. The skin friction coefficient C f and the local Nusselt number Nu x are defined as [43]: By employing Equation (7), one gets:

Stability Analysis
The temporal stability of the dual solutions as time evolves is studied. This analysis was first introduced by Merkin [44] and then followed by Weidman et al. [45] Firstly, consider the new variables as follows: Now, the unsteady form of Equations (2) and (3) are considered, while Equation (1) remains unchanged. On using (15), one obtains: subject to: Then, consider the following perturbation functions [45]: Here, Equation (19) is introduced to apply a small disturbance on the steady solution f = f 0 (η) and θ = θ 0 (η) of Equations (9)- (11). The functions F(η) and G(η) in Equation (19) are relatively small compared to f 0 (η) and θ 0 (η). The sign (positive or negative) of the eigenvalue γ determines the stability of the solutions. By employing (19), Equations (16)- (18) become: 1 Pr subject to: Without loss of generality we set F (0) = 1 [46] to get the eigenvalues γ in Equations (20) and (21). The stability of the solutions as time evolves is determined by examining the values of the smallest eigenvalue that was obtained. As time passes, there is an initial decay of disturbance if γ is positive (see Equation (19)), and thus the solution is stable and physically reliable in the long run. On the other hand, if γ is negative, there is an initial growth of disturbance, hence the solution is unstable.

Results and Discussion
Equations (9)-(11) were solved numerically by utilising the boundary value problem solver (bvp4c) in MATLAB software, which employs the 3-stage Lobatto IIIa formula [47]. This is a collocation formula and provides a continuous solution with fourth-order accuracy. The effectiveness of this solver ultimately counts on our ability to provide the algorithm with an initial guess for the solution. Moreover, the suitable value of the boundary layer thickness must be chosen depending on the values of the parameters applied. To solve this boundary value problem, it is necessary to first reduce the equations to a system of first-order ordinary differential equations. The effects of the physical parameters such as Al 2 O 3 (ϕ 1 ) and Cu (ϕ 2 ) nanoparticles volume fractions, the Prandtl number Pr, the radiation parameter R, and the mixed convection parameter λ on the flow behaviour are examined.
The values of the skin friction coefficient f (0) and the local Nusselt number −θ (0) for several values of Pr when R = 0, λ = 1, and ϕ 1 = ϕ 2 = 0 (regular fluid) are compared with published results of Ishak et al. [5], as presented in Table 3. It should be mentioned that Ishak et al. [5] solved their problem by the Keller-box method. Meanwhile, the boundary value problem solver (bvp4c) is employed in this study. It is found that the results are in excellent agreement. This gives confidence to the validity and accuracy of the numerical results for other values of parameters. Besides, the values of f (0) show a decreasing behaviour, while the values of −θ (0) increase for larger Pr. Additionally, Table 4 describes the values of Re 1/2 x C f and Re −1/2 x Nu x for Cu/water nanofluid when ϕ 1 = R = 0 and Pr = 6.2 with different values of λ and ϕ 2 . Here, we note that the values of both Re 1/2 x C f and Re −1/2 x Nu x increase with the increasing of λ and ϕ 2 . Besides, dual solutions are found for opposing (λ = −1) and assisting (λ = 1) flows, whereas the unique solution is obtained for λ = 0 (force convection flow). Furthermore, the values of Re 1/2 x C f for λ = 0 provided in the same table are compared with those of Bachok et al. [14], and the results are in excellent agreement, which thus gives confidence to the results for other values of λ. Table 3. Values of f (0) and −θ (0) for different values of Pr when ϕ 1 = ϕ 2 = 0 (regular fluid), R = 0, and λ = 1.
Pr.  Moreover, Table 5 shows the effect of λ, R and ϕ 2 on Re 1/2 x C f and Re −1/2 x Nu x when Pr = 6.2 for nanofluid (Cu/water) and hybrid nanofluid (Cu-Al 2 O 3 /water). For the first solutions, we found that the values of Re 1/2 x C f are accelerated with the increasing of λ and ϕ 2 ; however, they are decelerated with R. Besides, the values of Re −1/2 x Nu x enhance with increasing values of these parameters. The local Nusselt number Re −1/2 x Nu x enhance up to 32.37% for Cu-Al 2 O 3 /water (ϕ 1 = 0.1, ϕ 2 = 0.04) compared to the regular fluid (ϕ 1 = ϕ 2 = 0) when λ = −1, R = 0, and Pr = 6.2. Meanwhile, the values of Re −1/2 x Nu x are prominent for larger radiation (R = 1) with 49.63% enhancement compared to the non-radiant case (R = 0) when λ = −1, ϕ 1 = 0.1, ϕ 2 = 0.04, and Pr = 6.2. Moreover, the rise in λ from −1 to 1 contributes to the increment in the values of Re −1/2 x Nu x up to 8.66% when R = 1, ϕ 1 = 0.1, ϕ 2 = 0.04, and Pr = 6.2. The variations of the skin friction coefficient Re 1/2 x C f and the local Nusselt number Re −1/2 x Nu x against λ for several values of ϕ 2 and R are illustrated in Figures 2-5. The dual solutions of Equations (9)-(11) are possible for both assisting (λ > 0) and opposing (λ < 0) flows. The flow is accelerated for λ > 0 because there is a favourable pressure gradient induced by the buoyancy forces, which results in larger heat transfer and skin friction coefficients rather than the case of λ = 0 (non-buoyant case). We note that the separation of the boundary layer occurs when λ < 0. The dual solutions happen for λ > λ c and no solution for λ < λ c . The curve terminates at λ = λ c (critical value) and this point is known as the bifurcation point of the solutions. Separately, Figures 2 and 3 display the variations of Re 1/2 x C f and Re −1/2 x Nu x against λ for different values of ϕ 2 when Pr = 6.2 and ϕ 1 = 0.1 in the absence of R. It is observed that the values of Re 1/2 x C f and Re −1/2 x Nu x enhance with the rising of ϕ 2 . Moreover, it is noticed that the boundary layer separation is delayed with the added hybrid nanoparticles. The critical values are λ c = −4.6983, −5.1215, and −5.5404 for ϕ 2 = 0, 0.02 and 0.04, respectively. Apart from that, the variations of Re 1/2 x C f and Re −1/2 x Nu x with λ for different values of R when Pr = 6.2, ϕ 1 = 0.1, and ϕ 2 = 0.04 are illustrated in Figures 4 and 5. In the presence of R, we found that the skin friction coefficient Re 1/2 x C f decreases for λ < 0 but increases for λ > 0, whereas the local Nusselt number Re −1/2 x Nu x enhances for both cases. Besides, we notice that the domain of λ for the existence of the dual solutions decreases for larger values of R where the critical values of λ slightly increase. Note that the critical values λ c for R = 0, 1, and 2 are λ c = −5.5404, −4.7843, and −4.4030, respectively. It is observed in Figures 3 and 5, the second solutions of Re −1/2 x Nu x are boundless as λ → 0 − and as λ → 0 + .          The impact of ϕ 2 and R on the velocity f (η) and the temperature θ(η) profiles for the case of the opposing (λ = −1) and assisting (λ = 1) flows are presented in Figures 6-13. There exist dual solutions for f (η) and θ(η) which satisfy the infinity boundary conditions (11) asymptotically. The rising of ϕ 2 leads to an upsurge in the values of f (η) and θ(η) on the first solutions for both cases when Pr = 6.2, ϕ 1 = 0.1, and R = 0 as shown in Figures 6-9. Meanwhile, the velocity f (η) decreases when λ = −1 but increases when λ = 1 on the first solutions for larger values of R when Pr = 6.2, ϕ 1 = 0.1, and ϕ 2 = 0.04. The effect of R is to increase the temperature θ(η) inside the boundary layer for both cases as displayed in Figures 10-13. The radiation is dominant over conduction for larger values of R, causing a rise in the fluid temperature. It is also noticed that the solutions of the lower branch for the velocity have negative values ( f (η) < 0), which implies that the reverse flow occurs away from the wall, and these behaviors are displayed in Figures 6, 8, 10 and 12. The behaviors of θ(η) with different values of ϕ 2 and R for both cases are given in Figures 7, 9, 11 and 13. The overshoot of the temperature θ(η) near the wall is observed when λ = −1, and θ(η) < 0 when λ = 1 for the second solution.
The variations of γ against λ when Pr = 6.2, ϕ 1 = 0.1, ϕ 2 = 0.04 and R = 1 are described in Figure 14. For positive values of γ, it is noted that e − γτ → 0 as time evolves ( τ → ∞ ). In the meantime, for the negative value of γ, e − γτ → ∞ . These behaviors show that the first solution is stable and physically reliable, while the second solution is unstable in the long run. Mathematics 2021, 9, x FOR PEER REVIEW 10 of 16

Conclusions
The stagnation point flow towards a vertical plate in a hybrid nanofluid with thermal radiation was examined in the present paper. Findings revealed that dual solutions appeared for both assisting ( > 0) and opposing ( < 0) flows. The dual solutions were found for > and no solution for < , while the solutions bifurcated at = . In addition, the consequence of the copper nanoparticle volume fractions is to enhance the skin friction coefficient / and the local Nusselt number / for both cases. However, the values of / decreased for < 0, but increased for > 0, whereas the values of / were intensified for both cases in the presence of the radiation parameter . From these findings, the increments of the local Nusselt number / are observed in the range of 8.66% to 49.63% for the pertinent physical parameters considered. Besides, we noticed that the domain of the mixed convection parameter where the dual solutions are in existence decreased for larger . Further, the first solution of the velocity ′( ) and the temperature ( ) profiles enlarged with the increase of the copper nanoparticles volume fractions . Moreover, the effect of the radiation parameter is to increase the temperature ( ) inside the boundary layer for both cases. Lastly, it was discovered that between the two solutions, the solution with lower boundary layer thickness is stable and thus physically reliable in the long run.

Conclusions
The stagnation point flow towards a vertical plate in a hybrid nanofluid with thermal radiation was examined in the present paper. Findings revealed that dual solutions appeared for both assisting (λ > 0) and opposing (λ < 0) flows. The dual solutions were found for λ > λ c and no solution for λ < λ c , while the solutions bifurcated at λ = λ c . In addition, the consequence of the copper nanoparticle volume fractions ϕ 2 is to enhance the skin friction coefficient Re 1/2 x C f and the local Nusselt number Re −1/2 x Nu x for both cases. However, the values of Re 1/2 x C f decreased for λ < 0, but increased for λ > 0, whereas the values of Re −1/2 x Nu x were intensified for both cases in the presence of the radiation parameter R. From these findings, the increments of the local Nusselt number Re −1/2 x Nu x are observed in the range of 8.66% to 49.63% for the pertinent physical parameters considered. Besides, we noticed that the domain of the mixed convection parameter λ where the dual solutions are in existence decreased for larger R. Further, the first solution of the velocity f (η) and the temperature θ(η) profiles enlarged with the increase of the copper nanoparticles volume fractions ϕ 2 . Moreover, the effect of the radiation parameter R is to increase the temperature θ(η) inside the boundary layer for both cases. Lastly, it was discovered that between the two solutions, the solution with lower boundary layer thickness is stable and thus physically reliable in the long run. Acknowledgments: The authors would like to thank the anonymous reviewers for their constructive comments and suggestions. The financial supports received from the Universiti Kebangsaan Malaysia (Project Code: DIP-2020-001) and the Universiti Teknikal Malaysia Melaka are gratefully acknowledged.