Hybrid Nanofluid Flow over a Permeable Shrinking Sheet Embedded in a Porous Medium with Radiation and Slip Impacts

Abstract

The study of hybrid nanofluid and its thermophysical properties is emerging since the early of 2000s and the purpose of this paper is to investigate the flow of hybrid nanofluid over a permeable Darcy porous medium with slip, radiation and shrinking sheet. Here, the hybrid nanofluid consists of Cu/water as the base nanofluid and Al₂O₃–Cu/water works as the two distinct fluids. The governing ordinary differential equations (ODEs) obtained in this study are converted from a series of partial differential equations (PDEs) by the appropriate use of similarity transformation. Two methods of shooting and bvp4c function are applied to solve the involving physical parameters over the hybrid nanofluid flow. From this study, we conclude that the non-uniqueness of solutions exists through a range of the shrinking parameter, which produces the problem of finding a bigger solution than any other between the upper and lower branches. From the analysis, one can observe the increment of heat transfer rate in hybrid nanofluid versus the traditional nanofluid. The results obtained by the stability of solutions prove that the upper solution (first branch) is stable and the lower solution (second branch) is not stable.

1. Introduction

One of the most important industrial processes is heat transfer, carried out by heat exchangers in single and multiphase flow applications. Much interest and effort has created for experimental work in heat transfer due to the necessary need and solid demand for industrial applications that require the optimization and design of heat exchangers, despite the well-developed and built-in theoretical models that have existed since the 1970s. Many attempts have been made within these past years to enhance heat transfer rate, and one of the methods is by increasing the thermal conductivity. Choi [1] pioneered the first work of nanofluid and its capability in suspending nanoscale particles in the base fluid since they exhibit enhanced thermal conductivity and coefficient of convective heat transfer. Nanofluids have novel properties that make them a phenomenal development in many industrial applications, including microelectronics, hybrid-powered engines, domestic refrigerator, chiller, and even in high-functional military specialized gadgets, as explained by Saidur et al. [2]. One of the recent application in nanofluids was presented by Moghadasi et al. [3], who investigated the efficiency of synthesized nanosilica particles in reducing fines migration in hydraulic fracturing. They stated that the hydraulic fracturing process can be badly affected by fines migration, and they conducted an experiment by adding the nanosilica particles and soaking time. Their results have led to the clearer effluent fluid and less concentrations of clay particles in solution.

In these recent twenty years, many studies on nanofluids have been presented, and most of the reviews agree that one of the methods to change the pattern of flow is by considering the inclusion of nanoparticles into the base fluid. Motsumi and Makinde [4] performed a study on boundary layer flow over a permeable moving plate with nanofluids, viscous dissipation, and thermal radiation. Here, they considered Al₂O₃ and Cu as two distinct nanoparticles, and the results reveal that Al₂O₃ shows a higher velocity and thermal boundary layer than Cu. Sheikholeslami [5] studied the nanofluid flow and heat transfer over a cylinder with a uniform suction and described the increasing function of Nusselt number alongside nanoparticle volume fraction. Later, a study of magnetohydrodynamics (MHD) flow over a permeable stretching/shrinking sheet with nanofluid and suction/injection was published by Naramgari and Sulochana [6]. Based on their study, they indicated that the magnetic field parameter reduces the boundary layer flow, friction factor, and heat transfer rate on stretching surface. Other works that can be considered are found in [7,8,9,10,11,12,13,14,15].

The simplest way of defining a porous media is as a material that contains passages and filled with flowing fluid in liquid or gaseous forms. Examples of porosity are intergranular and intercrystalline, which are identified by their differences in molecular and cavern interstices. Hence, the potential of porous media has attracted much consideration in processing applications, as well as in academics works and publications. Ahmad and Pop [16] investigated the mixed convection flow through a vertical flat plate filled with nanofluids and porosity from a porous media, and they reported that two branches appeared, which are termed the lower and upper branches defined to the curves where the critical point of mixed convection parameter occurs. Sheikholeslami et al. [17] worked on nanofluid in porous media with magnetoyhydrodynamic transportation. In this study, they considered CuO-water as the nanoparticles in a porous cavity, and models of Darcy and Koo–Kleinstreuer–Li approach (KKL) were used to solve porous media and nanofluid, respectively. Bakar et al. [18] studied the mixed convection through a cylinder with nanofluid and thermal radiation saturated in a porous media, and they concluded that nanoparticles of alumina showed the highest rate in separating the boundary layer thickness, followed by titanium and copper. The studies of nanofluid in porous media are also successfully reported in [19,20,21,22].

Recently, researchers have gained much interest in hybrid nanofluid since numerous reports claim that the new hybrid nanoparticles may improve the heat transfer rate versus the classic nanofluid as well as minimize production cost, and these advantages can achieve a successful production for organizations, researchers, and academicians, as explained by Ghadikolaei et al. [23]. A hybrid nanofluid can be elaborated as a mixture of two or more different materials of nanometer sizes. Sundar et al. [24] classified hybrid nanofluids as a motivation in preparing the fluid flow to obtain further increment of heat transfer rate with augmented thermal conductivity of the involving nanofluids. In recent years, the authors of [25,26,27,28,29] conducted several other studies on hybrid nanofluid. Khashi’ie et al. [30] numerically studied the mixed convection of Cu–Al₂O₃/water in a non-Darcy porous medium with thermal dispersion. They found that Cu–Al₂O₃/water has greater heat transfer rate than nanofluid and regular fluid for some of the investigated parameters. In addition, a considerable amount of previous works on hybrid nanofluid over a porous medium have been successfully reported (e.g., [31,32,33,34,35]).

The numerical solution of Al₂O₃–Cu/water hybrid nanofluid along a permeable Darcy porous medium is conducted in this present work as the authors are inspired by the above-mentioned literature. We consider shrinking surface, slip factor, and radiation effect in this model. The main objective of this paper is to find the solutions to the current problem, which may benefit other researchers or academicians from the final outcomes.

2. Problem Formulation

In this paper, we contemplate a two-dimensional flow and heat transfer in a standard form over a permeable Darcy porous media with slip, hybrid nanofluid, radiation, and shrinking sheet, as illustrated in Figure 1. The coordinates of Cartesian on x and y are built-in along the sheet surface, and the surface is located at $y=0$. The x axis is chosen to be parallel in the direction of the surface motion, while the y axis is chosen to be perpendicular to the x axis.

Here, we review Copper (Cu) and Aluminium Oxide (Al₂O₃) as nano-sized particles and water as a base fluid. Table 1 lists the nanofluids and hybrid nanofluids thermophysical properties. We consider Cu and Al₂O₃ in this study as we follow the model introduced by [36], since these two nanoparticles are the most commonly used by many researchers in their experiment works and theoretical studies. It is noted that the basic thermophysical properties of nanofluid are extracted from the standard literature, and their properties of suspended nanoparticles versus fluid at 25° are listed in Table 2. We apply the Darcy equation in this model as it describes the fluid flow over a porous media, as suggested by Rajagopal [37]. Under the above assumptions, the continuity, momentum, and energy of nanoparticles equations based on Darcy flow model (see [38]) are as follows.

(1)${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=0,$

(2)$u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}-U_e{\displaystyle \frac{\partial U_e}{\partial x}}-{\nu }_{hnf}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+{\displaystyle \frac{{\nu }_{hnf}ϵ}{K}}(u-U_e)=0,$

(3)$u{\displaystyle \frac{\partial T}{\partial x}}+v{\displaystyle \frac{\partial T}{\partial y}}-{\displaystyle \frac{k_{hnf}}{{\left(\rho C_p\right)}_{hnf}}}{\displaystyle \frac{{\partial }^{2}T}{\partial y^2}}+{\displaystyle \frac{\partial q_r}{\partial y{\left(\rho C_p\right)}_{hnf}}}=0.$

Considering Rosseland’s approximation for radiation, as proposed by Rosseland [39] and Motsumi and Makinde [4], we have the radiative heat flux $q_r$ at

(4)$q_r=-{\displaystyle \frac{4\sigma }{3k\ast }}{\displaystyle \frac{\partial T^4}{\partial y}}.$

$T^4$ may be expressed as a temperature linear function and can be expanded using a truncated Taylor series since the difference in temperature is relatively small within the flow, for which we get $T^4\cong 4T_{\infty }^{3}T-3T_{\infty }^4$ by expanding $T^4$ and $T_{\infty }$. Hence, Equation (3) now can be reduced to

(5)$u{\displaystyle \frac{\partial T}{\partial x}}+v{\displaystyle \frac{\partial T}{\partial y}}={\displaystyle \frac{k_{hnf}}{{\left(\rho C_p\right)}_{hnf}}}{\displaystyle \frac{{\partial }^{2}T}{\partial y^2}}-{\displaystyle \frac{16\sigma T_{\infty }^3}{3k\ast {\left(\rho C_p\right)}_{hnf}}}{\displaystyle \frac{{\partial }^{2}T}{\partial y^2}}.$

The boundary conditions are now given by

(6)$\begin{array}{c}\hfill u=cx+L_1{\displaystyle \frac{\partial u}{\partial y}},v=v_w,T=T_w+D_1{\displaystyle \frac{\partial T}{\partial y}}\mathrm{at}y=0,\\ \hfill u\to U_e,T\to T_{\infty }\mathrm{as}y\to \infty .\end{array}$

The velocity components for the hybrid nanofluid along x and y axes are aligned with u and v, respectively; the hybrid nanofluid temperature is T; the external flow velocity is $U_e$ where $U_e=ax$; the porous media permeability is K; the dimensionless porosity of porous media is $ϵ$; the mean absorption rate of the nanofluid is $k\ast$; the constant number of Stefan–Boltzmann is $\sigma$; and $k_{hnf}$, ${\nu }_{hnf}$, ${\rho }_{hnf}$, and ${\left(\rho C_p\right)}_{hnf}$ are the hybrid nanofluids thermal conductivity, dynamic viscosity, density, and heat capacity, respectively. From Equation (6), the stretching/shrinking constant is c, the straining rate parameter is a, the suction or injection velocity constant is $v_w$, and $L_1$ and $D_1$ are the velocity and thermal slip factors, respectively.

Following Devi and Devi [36], the stream function and similarity transformations are introduced by

(7)$u=ax{f}^{\prime }\left(\eta \right)v=-\sqrt{a{\nu }_f}f\left(\eta \right),\theta \left(\eta \right)={\displaystyle \frac{T-T_{\infty }}{T_w-T_{\infty }}},\eta =y\sqrt{{\displaystyle \frac{a}{{\nu }_f}}}.$

By invoking the similarity variables in Equation (7) into Equations (2)–(5), we now have the new model of ODEs as follows

(8)${\displaystyle \frac{1}{A_1}}{f}^{‴}+f{f}^{″}-{\left(f^{\prime }\right)}^2-m_1{f}^{\prime }+m_1+1=0,$

(9)$(A_2+{\displaystyle \frac{4}{3}}R){\theta }^{″}+Prf{\theta }^{\prime }=0,$

subject to the boundary conditions at

(10)$\begin{array}{c}\hfill f\left(0\right)=S,f^{\prime }\left(0\right)=\alpha +\delta f^{″}\left(0\right),\theta \left(0\right)=1+\beta {\theta }^{\prime }\left(0\right),\\ \hfill f^{\prime }\left(\eta \right)\to 1,\theta \left(\eta \right)\to 0\mathrm{as}\eta \to \infty .\end{array}$

Here, the porous media permeability parameter is $m_1={\displaystyle \frac{{\nu }_{nf}ϵ}{Ka}}$, the radiation parameter is $R={\displaystyle \frac{4\sigma T_{\infty }^3}{k_{f}k\ast }}$, the local Prandtl number is $Pr={\displaystyle \frac{{\mu }_f{C}_{pf}}{k_f}}$, the suction parameter is $S={\displaystyle \frac{{\nu }_0}{\sqrt{c{\nu }_f}}}$, the shrinking parameter is $\alpha ={\displaystyle \frac{c}{a}}$, the slip velocity parameter is $\delta ={\displaystyle \frac{L_1{U}_{\infty }}{{\nu }_f}}$, and the slip thermal parameter is $\beta ={\displaystyle \frac{D_1{U}_{\infty }}{{\nu }_f}}$. The constants of $A_1$ and $A_2$ from Equations (8) and (9) are elaborated by

(11)$A_1={\displaystyle \frac{{\nu }_f}{{\nu }_{nf}}}={(1-{\varphi }_1)}^{2.5}{(1-{\varphi }_2)}^{2.5}\left\{{\varphi }_2\right({\displaystyle \frac{{\rho }_{s1}}{{\rho }_f}})+(1-{\varphi }_2)[(1-{\varphi }_1)+{\varphi }_1\left({\displaystyle \frac{{\rho }_{s1}}{{\rho }_f}}\right)\left]\right\},$

(12)$A_2={\displaystyle \frac{k_{hnf}/k_f}{\{(1-{\varphi }_2)[(1-{\varphi }_1)+{\varphi }_1({\displaystyle \frac{{\left(\rho C_p\right)}_{s1}}{{\left(\rho C_p\right)}_{s2}}}\left)\right]+{\varphi }_2\left({\displaystyle \frac{{\left(\rho C_p\right)}_{s2}}{{\left(\rho C_p\right)}_{s1}}}\right)\}}}.$

The current study requires physical quantities of interest which are skin friction coefficient $C_f$ and the local Nusselt number $N{u}_x$. Hence, the responding $C_f$ and $N{u}_x$ are

(13)$C_f={\displaystyle \frac{{\tau }_w}{{\rho }_f{U}^2}},N{u}_x={\displaystyle \frac{x{q}_w}{k_f(T_w-T_{\infty })}},$

and, by simplifying Equations (7) and (13), we have

(14)$C_{fx}R{e}_x^{1/2}={\displaystyle \frac{f^{″}\left(0\right)}{{(1-{\varphi }_1)}^{2.5}{(1-{\varphi }_2)}^{2.5}}},N{u}_{x}R{e}_x^{-1/2}=-({\displaystyle \frac{k_{hnf}}{k_f}}+{\displaystyle \frac{4}{3}}R){\theta }^{\prime }\left(0\right).$

Here, the local Reynolds number is represented by $R{e}_x={\displaystyle \frac{U_{w}x}{{\nu }_f}}$.

3. Numerical Soluion

The system of ODEs in Equations (8) and (9) subjected to the boundary conditions in Equation (10) were numerically solved using the method of shooting technique via Maple and bvp4c function implemented in MATLAB (see [40]), with various numbers for different parameters. In numerical analysis, the shooting method is a technique for reducing a boundary value problem into a set of initial value problems in order to solve the problem. The method can be successfully achieved by shooting the trajectories in different directions until the desired boundary value has been found. Another way, bvp4c describes a finite difference code that employs the three-stage Lobatto Illa formula, as highlighted by Zainal et al. [41]. The bvp4c function is a collocation formula which provides the polynomial at a C⁻¹-continuous solution that is fourth-order accurate in the specific interval. Hence, the variable ${\eta }_{max}$ is acquired by applying the boundary conditions of the field at the finite value for the similarity variable $\eta$. Thus, we set ${\eta }_{max}=9$ in our analysis to fulfill the far field boundary conditions as in Equation (10) asymptotically.

Due to the eligibility and accuracy of our numerical result, a comparison is made among the present skin friction coefficient with those of Wang [42] and Bhattacharyya et al. [43], as shown in Table 3. Here, the parameters of $m_1$, S, R, $\delta$, and $\beta$ are absent, while Pr is standardized at 0.7. In the comparison, we observed a good agreement between the present and previous works.

Figure 2, Figure 3 and Figure 4 illustrate the dual solutions obtained from the skin friction coefficient $C_{fx}R{e}_x^{1/2}$ and local Nusselt number $N{u}_{x}R{e}_x^{-1/2}$ with various values of shrinking parameter $\alpha$, suction parameter S, velocity slip parameter $\delta$, and radiation parameter R. Figure 2 shows that the impact of $\alpha =1.0$ resulting in $C_{fx}R{e}_x^{1/2}=0$. This can be explained by the fact that no friction exists at the fluid–solid interface when the fluid and solid boundaries move at the same velocity. At the same time, a negative value emerges when $\alpha >1$, which indicates that a drag force is applied by the fluid along the boundary of solid, and vice versa. Figure 2 and Figure 3 show the increment of critical point in $C_{fx}R{e}_x^{1/2}$ and $N{u}_{x}R{e}_x^{-1/2}$ when the values of S and $\delta$ increase. The main reason for all these physical behaviors can be explained by the combination effects between shrinking sheet strength and porosity at the surface.

Figure 5 illustrates the numbers of volume particle parameter ${\varphi }_1$ for Al₂O₃ against velocity profiles $f^{\prime }\left(\eta \right)$ and temperature profiles $\theta \left(\eta \right)$. In these figures, we depict that the upper solution in $f^{\prime }\left(\eta \right)$ decreases, while the rest shows promising positive pattern along the flow. These behaviors of increase and decrease can be explained by a contribution of the flow and the conditions of thermal and dispersive elements properties that maximize the heat transfer.

We present the velocity profiles $f^{\prime }\left(\eta \right)$ and temperature profiles $\theta \left(\eta \right)$ for several values of porous media permeability parameter $m_1$ and suction parameter S in Figure 6 and Figure 7, respectively, where the behavior of $f^{\prime }\left(\eta \right)$ shows an increment in upper solution and a decrement in lower solution versus the increased number of $m_1$ and S. Meanwhile, the flow behavior of $\theta \left(\eta \right)$ depicts a reverse pattern where the solution shows a decrement.

To analyze the influence of the nanoparticles addition on the fields of thermal and dynamic flow, Figure 8 displays the streamlines for several numbers of volume particle parameter $\varphi$. By increasing the number of ${\varphi }_1$ for Al₂O₃/water and ${\varphi }_2$ for Cu/water, we note that the strength of flow increases, as can be seen from the pattern of the streamlines by alerting that the increase of Cu/water nanoparticle number has a higher heat transfer rate as compared to Al₂O₃/water nanofluid.

4. Stability Analysis

In this respect, the dual nature of the solutions is observed from our previous analysis, and, hence, it is necessary to perform a stability analysis in order to identify the stability of each solutions, as suggested by Merkin [44] and Merrill et al. [45]. Here, we consider our model of momentum and energy in an unsteady state and we have

(15)${\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}-U_e{\displaystyle \frac{\partial U_e}{\partial x}}-{\nu }_{hnf}{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+{\displaystyle \frac{{\nu }_{hnf}ϵ}{K}}(u-U_e)=0,$

(16)${\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial T}{\partial x}}+v{\displaystyle \frac{\partial T}{\partial y}}-{\displaystyle \frac{k_{hnf}}{{\left(\rho C_p\right)}_{hnf}}}{\displaystyle \frac{{\partial }^{2}T}{\partial y^2}}+{\displaystyle \frac{\partial q_r}{\partial y{\left(\rho C_p\right)}_{hnf}}}=0,$

where t represents time. Our boundary conditions now changes to

(17)$\begin{array}{c}\hfill u=v=0,T=T_{\infty }\mathrm{when}t<0,\\ \hfill u=cx+L_1{\displaystyle \frac{\partial u}{\partial y}},v=v_w,T=T_w+D_1{\displaystyle \frac{\partial T}{\partial y}}\mathrm{at}y=0,\\ \hfill u\to U_e,T\to T_{\infty }\mathrm{as}y\to \infty \mathrm{when}t\ge 0.\end{array}$

We now introduce a new dimensionless variable $\tau$ in regards of t, where we have

(18)$u=ax{f}^{\prime }\left(\eta \right)v=-\sqrt{a{\nu }_f}f\left(\eta \right),\theta \left(\eta \right)={\displaystyle \frac{T-T_{\infty }}{T_w-T_{\infty }}},\eta =\sqrt{{\displaystyle \frac{a}{{\nu }_f}}}y,\tau =at,$

so that Equations (15) and (16) can be formed into

(19)${\displaystyle \frac{1}{A_1}}{\displaystyle \frac{{\partial }^{3}f}{\partial {\eta }^3}}-{\left({\displaystyle \frac{\partial f}{\partial \eta }}\right)}^2+f{\displaystyle \frac{{\partial }^{2}f}{\partial {\eta }^2}}-m_1{\displaystyle \frac{\partial f}{\partial \eta }}+m_1+1-{\displaystyle \frac{{\partial }^{2}f}{\partial \eta \partial \tau }}=0,$

(20)$(A_2+{\displaystyle \frac{4}{3}}R){\displaystyle \frac{{\partial }^2\theta }{\partial {\eta }^2}}+Prf{\displaystyle \frac{\partial \theta }{\partial \eta }}-{\displaystyle \frac{\partial \theta }{\partial \tau }}=0,$

with respect to

(21)$\begin{array}{c}\hfill f(0,\tau )=S,{\displaystyle \frac{\partial f}{\partial \eta }}(0,\tau )=\alpha +\delta {\displaystyle \frac{{\partial }^{2}f}{\partial {\eta }^2}},\theta (0,\tau )=1+\beta {\displaystyle \frac{\partial \theta }{\partial \eta }}\mathrm{at}\eta =0,\\ \hfill {\displaystyle \frac{\partial f}{\partial \eta }}(\eta ,\tau )\to 1,\theta (\eta ,\tau )\to 0\mathrm{as}\eta \to \infty .\end{array}$

In regards of our dual solutions, one can adopt the analyses suggested by Merkin [46] and Weidman et al. [47], which are as follows

(22)$\begin{array}{c}\hfill f(\eta ,\tau )=f_0\left(\eta \right)+e^{-\gamma \tau }{F}_0(\eta ,\tau ),\\ \hfill \theta (\eta ,\tau )={\theta }_0\left(\eta \right)+e^{-\gamma \tau }{G}_0(\eta ,\tau ).\end{array}$

$\gamma$ in Equation (22) is a parameter of unknown eigenvalue, while $F_0\left(\eta \right)$ and $G_0\left(\eta \right)$ are the small relatives of $f_0\left(\eta \right)$ and ${\theta }_0\left(\eta \right)$, respectively. Here, $\gamma$ is infamous for the decay or growth of a disturbance, where the smallest $\gamma$ in positive number represents the continuous decaying of disturbances, in which we can finalize the solution to be in a stable state, and vice versa. To test our numerical procedure, we simplify Equations (19), (20), and (22) as

(23)${\displaystyle \frac{1}{A_1}}{F}_0^{‴}+f_0{F}_0^{″}-2f_0^{\prime }{F}_0^{\prime }+f_0^{″}{F}_0-m_1{F}_0^{\prime }+\gamma F_0^{\prime }=0,$

(24)$(A_2+{\displaystyle \frac{4}{3}}R)G_0^{″}+Pr{f}_0{G}_0^{\prime }+Pr{\theta }_0^{\prime }{F}_0+\gamma G_0^{\prime }=0.$

Our boundary conditions now can simplify to

(25)$\begin{array}{c}\hfill F_0\left(0\right)=0,F_0^{\prime }\left(0\right)=0,G_0\left(0\right)=0,\\ \hfill F_0^{\prime }\left(\eta \right)\to 0,G_0\left(\eta \right)\to 0\mathrm{as}\eta \to \infty .\end{array}$

In regards of Equation (25), Harris et al. [48] suggested relaxing the boundary condition on $F_0\left(\eta \right)\to 0$ and $G_0\left(\eta \right)\to 0$ for a fixed value of $\gamma$ in order to determine the range of possible eigenvalues. Here, we can solve the problem with the new boundary conditions at $F_0^{″}\left(\eta \right)\to 1$ as $\eta \to \infty$. Due to this formulation, we could analyze the stability of our dual solutions via bvp4c function in MATLAB software. The value of the smallest eigenvalues $\gamma$ against various ${\varphi }_1$ and ${\varphi }_2$ numbers are presented in Table 4 and Table 5, respectively. In both tables, it is noticed that a series of positive numbers appears throughout the upper solution (first branch), while a series of negative numbers is observed throughout the lower solution (second branch). Hence, a conclusion can be finalized that the first solution is stable and realizable, and vice versa.

5. Conclusions

A numerical investigation of a steady, two-dimensional hybrid nanofluid over a Darcy porous media past a permeable shrinking sheet with radiation is studied in this present work. A new type of Al₂O₃–Cu/water is employed in this study as a model of hybrid nanofluid. From our observation, we conclude that the skin friction coefficient $C_{fx}R{e}_x^{1/2}$ is expanding when we increase the number of involving parameters and most of the parameters used in this investigation show an increasing pattern on boundary layer flow, in either upper or lower solution. Moreover, two branches of solutions are found to exist within a range of negative numbers in shrinking parameter $\alpha$. Due to this, the most stable solution between these two is identified via a work of stability analysis. It is then concluded that the first branch (upper solution) is stable and physically realizable, while the second branch (lower solution) is unstable.

Author Contributions

Investigation, S.A.B.; Methodology, S.A.B.; Supervision, N.M.A.; Validation, N.S.K.; Writing—original draft, S.A.B.; Writing—review and editing, N.M.A., N.S.K. and N.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors gratefully acknowledge the Ministry of Higher Education Malaysia for financial support through the Fundamental Research Grant Scheme (KPT FRGS/1/ 2019/STG06/UPM/02/3, Vot 5540309).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

Al2O3	Aluminium Oxide
bvp4c	Boundary value problem – fourth-order method
Cu	Copper
CuO	Copper Oxide
KKL	Koo–Kleinstreuer–Li approach
MHD	Magnetohydrodynamics
ODE	Ordinary differential equations
PDE	Partial differential equations

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The system of coordinate and model that moving in a shrinking sheet.

Figure 2. CfxRex1/2 and NuxRex−1/2 for Al2O3–Cu/water with S and α.

Figure 2. CfxRex1/2 and NuxRex−1/2 for Al2O3–Cu/water with S and α.

Figure 3. CfxRex1/2 and NuxRex−1/2 for Al2O3–Cu/water with δ and α.

Figure 4. Local Nusselt number NuxRex−1/2 for Al2O3–Cu/water with R and α.

Figure 5. The flow of f′(η) and θ(η) against ϕ1.

Figure 6. The flow of f′(η) and θ(η) against m1.

Figure 7. The flow of f′(η) and θ(η) against S.

Figure 8. Streamlines for Cu and Al2O3/water.

Table 1

Properties of nanofluids and hybrid nanofluids.

Properties	Nanofluid	Hybrid Nanofluid
Density	${\rho }_{nf}=(1-\varphi ){\rho }_f+\varphi {\rho }_s$	${\rho }_{hnf}=(1-{\varphi }_2)[(1-{\varphi }_1){\rho }_f+{\varphi }_1{\rho }_{s1}]+{\varphi }_2{\rho }_{s2}$
Heat capacity	${\left(\rho C_p\right)}_{nf}=(1-\varphi ){\left(\rho C_p\right)}_f+\varphi {\left(\rho C_p\right)}_s$	${\left(\rho C_p\right)}_{hnf}=(1-{\varphi }_2)[(1-{\varphi }_1){\left(\rho C_p\right)}_f+$
		${\varphi }_1{\left(\rho C_p\right)}_{s1}]+{\varphi }_2{\left(\rho C_p\right)}_{s2}$
Dynamic viscosity	${\nu }_{nf}={\displaystyle \frac{{\nu }_f}{{(1-\varphi )}^{2.5}}}$	${\nu }_{hnf}={\displaystyle \frac{{\nu }_f}{{(1-{\varphi }_1)}^{2.5}{(1-{\varphi }_2)}^{2.5}}}$
Thermal conductivity	${\displaystyle \frac{k_f}{k_{nf}}}={\displaystyle \frac{k_s+2k_f+\varphi (k_f-k_s)}{k_s+2k_f-2\varphi (k_f-k_s)}}$	${\displaystyle \frac{k_{nf}}{k_{hnf}}}={\displaystyle \frac{k_{s2}+2k_{nf}+{\varphi }_2(k_{nf}-k_{s2})}{k_{s2}+2k_{nf}-2{\varphi }_2(k_{nf}-k_{s2})}}$
		where
		${\displaystyle \frac{k_f}{k_{nf}}}={\displaystyle \frac{k_{s1}+2k_f+{\varphi }_1(k_f-k_{s1})}{k_{s1}+2k_f-2{\varphi }_1(k_f-k_{s1})}}$

Table 2

Fluid and nanoparticles thermophysical characteristics (see [36]).

Physical Characteristics	Water (f)	Al2O3 (s1)	Cu (s2)
Density, $\rho$ (kg/m${}^3$)	997.0	3970	8933
Thermal expansion, $\beta$ (K${}^{-1}$)	$21\times {10}^{-5}$	$0.85\times {10}^{-5}$	$1.67\times {10}^{-5}$
Thermal conductivity, k (W/m K)	0.6071	40	400
Thermal capacity, $C_p$ (J/kg K)	4180	765	385

Table 3

Comparison of $C_{fx}R{e}_x^{1/2}$ for the present study and those by Wang [42] and Bhattacharyya et al. [43].

$\mathit{\alpha }$	${\mathit{C}}_{\mathit{f}\mathit{x}}{\mathit{R}\mathit{e}}_{\mathit{x}}^{1/2}$
	Present Study	Bhattacharyya et al. [43]	Wang [42]
−0.50	1.49566	1.49655	1.49567
−0.625	1.52071	1.50715	–
−0.75	1.48929	1.48929	1.48930
−1.00	1.32882(0)	1.32881(0)	1.32282(0)
−1.15	1.08223(0.11670)	1.08223(0.11670)	1.08223(0.11670)
−1.20	0.93247(0.23364)	0.93247(0.23364)	–

Table 4

Smallest eigenvalue numbers of $\gamma$ against ${\varphi }_1$.

${\mathit{\varphi }}_1$	$\mathit{\alpha }$	Upper Solution	Lower Solution
0.2	−1.6	0.87412	−0.56313
	−1.7	0.90673	−0.56422
	−1.8	0.94051	−0.56540
0.3	−1.6	0.91186	−0.59207
	−1.7	0.93004	−0.61358
	−1.8	0.96728	−0.63776

Table 5

Smallest eigenvalue numbers of $\gamma$ against ${\varphi }_1$.

${\mathit{\varphi }}_2$	$\mathit{\alpha }$	Upper Solution	Lower Solution
0.1	−1.6	0.51620	−0.47993
	−1.7	0.55783	−0.48015
	−1.8	0.59022	−0.48154
0.2	−1.6	0.73326	−0.52197
	−1.7	0.74748	−0.52244
	−1.8	0.81905	−0.52293