Hybrid Nanofluid Flow Past a Permeable Moving Thin Needle

Abstract

The problem of a steady flow and heat transfer past a permeable moving thin needle in a hybrid nanofluid is examined in this study. Here, we consider copper (Cu) and alumina (Al₂O₃) as hybrid nanoparticles, and water as a base fluid. In addition, the effects of thermophoresis and Brownian motion are taken into consideration. A similarity transformation is used to obtain similarity equations, which are then solved numerically using the boundary value problem solver, bvp4c available in Matlab software (Matlab_R2014b, MathWorks, Singapore). It is shown that heat transfer rate is higher in the presence of hybrid nanoparticles. It is discovered that the non-uniqueness of the solutions is observed for a certain range of the moving parameter $\lambda$. We also observed that the bifurcation of the solutions occurs in the region of $\lambda <0$, i.e., when the needle moved toward the origin. Furthermore, we found that the skin friction coefficient and the heat transfer rate at the surface are higher for smaller needle sizes. A reduction in the temperature and nanoparticle concentration was observed with the increasing of the thermophoresis parameter. It was also found that the increase of the Brownian motion parameter leads to an increase in the nanoparticle concentration. Temporal stability analysis shows that only one of the solutions was stable and physically reliable as time evolved.

1. Introduction

The development of advanced heat transfer fluids has received considerable coverage from the researchers and scientists over the last few years. Regular fluids (ethylene glycol, oil, water) are commonly used in the industrial and engineering applications. However, the heat transfer rate of these fluids is limited due to weak thermal conductivity. Therefore, to resolve its deficiency, a single form of nanosized particles is applied to the above-mentioned fluids and is called ‘nanofluid’. This term was introduced by Choi and Eastman [1] for the first time in 1995. To name a few, the advantages of utilizing nanofluids filled in a rectangular enclosure have been examined by Khanafer et al. [2] and Oztop and Abu-Nada [3].

Nevertheless, ‘hybrid nanofluid’ was developed to upgrade a regular nanofluid’s thermal properties. It seems that Turcu et al. [4] and Jana et al. [5] are among the earliest researchers who considered the hybrid nanocomposite particles in their experimental studies. Hybrid nanofluid is an advanced fluid that incorporates more than one nanoparticle which has the capacity of raising the heat transfer rate because of the synergistic effects [6]. In addition, the desired heat transfer can be accomplished by combining or hybridizing the suitable nanoparticles [7].

Over the past few years, the boundary layer flow past a stretched or shrunk surface in a hybrid nanofluid has been extensively investigated. The research in this area was rapidly established due to its important applications in the industrial processes, for example, in paper production, extraction of polymer, artificial fiber and glass blowing. The investigations of the hybrid nanofluid flow over a stretched surface considering 2D and 3D flows were done by Devi and Devi [8,9]. In these studies, they discovered that the larger nanoparticle volume fractions contributed to the enhancement of the heat transfer rate. Furthermore, a new hybrid nanofluid’s thermophysical model was introduced in their studies. The validity of this new thermophysical model was achieved by comparing the results with the experimental data from Suresh et al. [10]. Moreover, Hayat and Nadeem [11] studied the problem of a hybrid nanofluid composed of Ag–CuO/water for the three-dimensional rotating flow. In addition, the dual nature of the flow past a stretched and shrunk surface in a hybrid nanofluid with temporal stability analysis was reported by Waini et al. [12]. They discovered that one of the solutions is unstable in the long run, while the other is stable and thus physically reliable. After that, various flows and heat transfer of a hybrid nanofluid were extended to the different aspects by Waini et al. [13,14,15,16,17,18,19,20].

The flow over a thin needle has become a topic of interest for the researchers because of its importance in industrial applications such as the hot wire anemometer, electronic devices and geothermal power generation. The thin needle is viewed as a body of revolution where its thickness is smaller compared to the boundary layer thickness. Historically, Lee [21] initially studied the viscous fluid flow past a thin needle. Then, this work was explored by the researchers [22,23,24,25,26,27,28,29] by considering various physical aspects. Furthermore, in the work of Ishak et al. [30], the dual solutions were obtained when the movement of the needle is in the opposite direction to the free stream. Then, a similar problem was studied by Waini et al. [31] but with the heat flux surface temperature. Apart from that, the effect of the nanoparticle on the flow over a thin needle was initiated by Grosan and Pop [32] and then extended by Soid et al. [33] to a moving needle in a nanofluid. Furthermore, the effect of the magnetic field on the flow past a moving vertical thin needle in a nanofluid was reported by Salleh et al. [34].

This study aims to examine the flow past a permeable moving thin needle in a hybrid nanofluid by employing Tiwari and Das [35] and Buongiorno [36] nanofluid models. Here, we consider water as the base fluid, while copper (Cu) and alumina (Al₂O₃) as the hybrid nanoparticles. Further, the effect of thermophoresis and Brownian motion are also considered. The results are obtained for several physical parameters and presented graphically and through tables. Also, the comparison results for limiting cases are done with previously published data.

2. Mathematical Formulation

The flow configuration past a permeable moving thin needle in a hybrid nanofluid is demonstrated in Figure 1. Here, $x$ and $r$ are the cylindrical coordinates with $x$- and $r$- being the axial and radial coordinates, respectively. The needle moves with a constant speed $U_{\infty }$ in a moving fluid with the same speed $U_{\infty }.$ Furthermore, the surface temperature $T_w$ and the ambient temperature $T_{\infty }$ are constant such that $T_w>T_{\infty }$, whereas $C_{\infty }$ is the ambient nanoparticle concentration with the normal flux of nanoparticles is zero at the wall. The shape of the nanoparticle is spherical, and its size is uniform, while the agglomeration is disregarded since the hybrid nanofluid is formed as a stable composite. After employing the usual boundary layer approximations, the governing equations of the hybrid nanofluid are given by (see Soid et al. [33]; Kuznetsov and Nield [37]):

$$\frac{\partial }{\partial x}\left(ru\right)+\frac{\partial }{\partial r}\left(rv\right)=0$$

(1)

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial r}=\frac{{\mu }_{hnf}}{{\rho }_{hnf}}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)$$

(2)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial r}=\frac{1}{r}\frac{k_{hnf}}{{\left(\rho C_p\right)}_{hnf}}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)+\sigma \left[D_B\frac{\partial C}{\partial r}\frac{\partial T}{\partial r}+\frac{D_T}{T_{\infty }}{\left(\frac{\partial T}{\partial r}\right)}^2\right]$$

(3)

$$u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial r}=\frac{D_B}{r}\frac{\partial }{\partial r}\left(r\frac{\partial C}{\partial r}\right)+\frac{D_T}{T_{\infty }}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)$$

(4)

subject to:

$$\begin{array}{c}u=U_{\infty }\lambda , v=v_w, T=T_w, D_B\frac{\partial C}{\partial r}+\frac{D_T}{T_{\infty }}\frac{\partial T}{\partial r}=0 \mathrm{at} r=R\left(x\right)\\ u\to U_{\infty }, T\to T_{\infty }, C\to C_{\infty } \mathrm{as} r\to \infty \end{array}$$

(5)

where the surface of the thin needle is described by $R\left(x\right)$, $u$ and $v$ are the velocity components in the axial $x$- and radial $r$- directions, $D_T$ and $D_B$ are thermophoretic diffusion and Brownian diffusion coefficients, respectively. Additionally, $\sigma , T, C$ and $v_w$ represent the effective heat capacity ratio, hybrid nanofluid’s temperature, nanoparticle concentration and mass flux velocity.

Further, ${\left(\rho C_p\right)}_{hnf},k_{hnf},{\mu }_{hnf},$ and ${\rho }_{hnf}$ characterize the heat capacity, thermal conductivity, dynamic viscosity and density of the hybrid nanofluid, respectively, where their thermophysical properties are defined in

Table 1. Thermophysical properties of nanofluid and hybrid nanofluid (see [3,8,15]).

Thermophysical Properties	Nanofluid	Hybrid Nanofluid
Density	${\rho }_{nf}=\left(1-{\phi }_1\right){\rho }_f+{\phi }_1{\rho }_{n1}$	${\rho }_{hnf}=\left(1-{\phi }_2\right)\left[\left(1-{\phi }_1\right){\rho }_f+{\phi }_1{\rho }_{n1}\right]+{\phi }_2{\rho }_{n2}$
Heat capacity	${(\rho C_p)}_{nf}=\left(1-{\phi }_1\right) {(\rho C_p)}_f+{\phi }_1 {(\rho C_p)}_{n1}$	${\left(\rho C_p\right)}_{hnf}=\left(1-{\phi }_2\right)\left[\left(1-{\phi }_1\right){\left(\rho C_p\right)}_f+{\phi }_1{(\rho C_p)}_{n1}\right]+{\phi }_2{\left(\rho C_p\right)}_{n2}$
Dynamic viscosity	${\mu }_{nf}=\frac{{\mu }_f}{{(1-{\phi }_1)}^{2.5}}$	${\mu }_{hnf}=\frac{{\mu }_f}{{\left(1-{\phi }_1\right)}^{2.5} {\left(1-{\phi }_2\right)}^{2.5}}$
Thermal conductivity	$\frac{k_{nf}}{k_f}=\frac{k_{n1}+2k_f-2{\phi }_1\left(k_f-k_{n1}\right)}{k_{n1}+2k_f+{\phi }_1\left(k_f-k_{n1}\right)}$	$\frac{k_{hnf}}{k_{nf}}=\frac{k_{n2}+2k_{nf}-2{\phi }_2\left(k_{nf}-k_{n2}\right)}{k_{n2}+2k_{nf}+{\phi }_2\left(k_{nf}-k_{n2}\right)}$ where $\frac{k_{nf}}{k_f}=\frac{k_{n1}+2k_f-2{\phi }_1\left(k_f-k_{n1}\right)}{k_{n1}+2k_f+{\phi }_1\left(k_f-k_{n1}\right)}$

. Meanwhile, the physical properties of Cu, Al₂O₃ and water are given in

Table 2. Thermophysical properties of nanoparticles and water (see [3,15]).

Thermophysical Properties	Al2O3	Cu	Water
$\rho \left(kg/m^3\right)$	3970	8933	997.1
$C_p \left(J/kgK\right)$	765	385	4179
$k \left(W/mK\right)$	40	400	0.613
Prandtl number, $\mathrm{Pr}$			6.2

. Here, $k,\rho ,\mu ,\left(\rho C_p\right),$ and $C_p$ represent the thermal conductivity, density, dynamic viscosity, heat capacity and specific heat at constant pressure, respectively. Meanwhile, the nanoparticle volume fractions of Al₂O₃ and Cu are symbolized by ${\phi }_1$ and ${\phi }_2$, while their solid components are indicated by the subscripts $n1$ and $n2$, respectively. In addition, the hybrid nanofluid, nanofluid and fluid, are indicated by the subscripts $hnf,nf$ and $f$, respectively.

An appropriate transformation is introduced as follows (see Soid et al. [33] and Kuznetsov and Nield [37]):

$$\psi ={\nu }_{f}x f\left(\eta \right), \theta \left(\eta \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }}, \varphi \left(\eta \right)=\frac{C-C_{\infty }}{C_{\infty }}, \eta =\frac{U_{\infty } r^2}{{\nu }_f x}$$

(6)

where $\psi$ denotes the stream function with $u=\left(1/r\right)\left(\partial \psi /\partial r\right)$ and $v=-\left(1/r\right)\left(\partial \psi /\partial x\right)$ so that Equation (1) is identically satisfied. Thus, we have:

$$u=2U_{\infty }{f}^{\prime }\left(\eta \right), v=-\frac{{\nu }_f}{r}\left(f\left(\eta \right)-\eta f^{\prime }\left(\eta \right)\right)$$

(7)

and,

$$v_w=-\frac{{\nu }_f}{r}\frac{\lambda }{2}c$$

(8)

where primes denote differentiation with respect to $\eta$ and ${\nu }_f$ is the fluid’s kinematic viscosity. Setting $\eta =c$ (refer to the wall of the needle) in Equation (6), the needle’s surface can be defined as:

$$R\left(x\right)={\left(\frac{{\nu }_{f}c x}{U_{\infty }}\right)}^{1/2}$$

(9)

Using (6), Equations (2) to (4) become:

$$2\frac{{\mu }_{hnf}/{\mu }_f}{{\rho }_{hnf}/{\rho }_f}{\left(\eta f^{″}\right)}^{\prime }+f{f}^{″}=0$$

(10)

$$\frac{2}{\mathrm{Pr}}\frac{k_{hnf}/k_f}{{(\rho C_p)}_{hnf}/{(\rho C_p)}_f}{\left(\eta {\theta }^{\prime }\right)}^{\prime }+f{\theta }^{\prime }+2\eta \left(Nb{\varphi }^{\prime }{\theta }^{\prime }+Nt{{\theta }^{\prime }}^2\right)=0$$

(11)

$$2{\left(\eta {\varphi }^{\prime }\right)}^{\prime }+Scf{\varphi }^{\prime }+2\frac{Nt}{Nb}{\left(\eta {\theta }^{\prime }\right)}^{\prime }=0$$

(12)

subject to:

$$\begin{array}{c}f\left(c\right)=\lambda c, f^{\prime }\left(c\right)=\frac{\lambda }{2}, \theta \left(c\right)=1, Nb{\varphi }^{\prime }\left(c\right)+Nt{\theta }^{\prime }\left(c\right)=0\\ f^{\prime }\left(\eta \right)\to \frac{1}{2}, \theta \left(\eta \right)\to 0, \varphi \left(\eta \right)\to 0 \mathrm{as} \eta \to \infty \end{array}$$

(13)

Here, $\mathrm{Pr}$ represents the Prandtl number, $Sc$ is the Schmidt number, $Nt$ and $Nb$ denote the thermophoresis and Brownian motion parameters, respectively and $\lambda$ is the constant moving parameter. The needle moves away from the origin if $\lambda >0$ and moves toward the origin if $\lambda <0$. All these parameters are defined as:

$$\mathrm{Pr}=\frac{{\left(\mu C_p\right)}_f}{k_f}, Sc=\frac{{\nu }_f}{D_B}, Nt=\frac{\sigma D_T \left(T_w-T_{\infty }\right)}{{\nu }_f T_{\infty }}, Nb=\frac{\sigma D_B C_{\infty }}{{\nu }_f}$$

(14)

The physical quantities of interest are the skin friction coefficient $C_f$, the local Nusselt number $N{u}_x$ and the local Sherwood number $S{h}_x$ which are defined as

$$C_f=\frac{{\tau }_w}{{\rho }_f U_{\infty }^2}, N{u}_x=\frac{x q_w}{k_f \left(T_w-T_{\infty }\right)}, S{h}_x=\frac{x q_m}{D_B C_{\infty }},$$

(15)

where the shear stress along the surface of the needle ${\tau }_w$, the heat flux $q_w$ and the mass flux from the surface of the needle $q_m$ are given by

$${\tau }_w={\mu }_{hnf}{\left(\frac{\partial u}{\partial r}\right)}_{r=R\left(x\right)}, q_w=-k_{hnf}{\left(\frac{\partial T}{\partial r}\right)}_{r=R\left(x\right)}, q_m=-D_B{\left(\frac{\partial C}{\partial r}\right)}_{r=R\left(x\right)}$$

(16)

Using (6), (15) and (16), we get

$$\begin{array}{c}R{e}_x^{1/2}{C}_f=4c^{1/2}\frac{{\mu }_{hnf}}{{\mu }_f}{f}^{″}\left(0\right), R{e}_x^{-1/2}N{u}_x=-2c^{\frac{1}{2}}\frac{k_{hnf}}{k_f}{\theta }^{\prime }\left(c\right)\\ R{e}_x^{-1/2}S{h}_x=-2c^{1/2}{\varphi }^{\prime }\left(c\right)=2c^{1/2}\frac{Nt}{Nb}\theta \prime \left(c\right),\end{array}$$

(17)

where $R{e}_x=U_{\infty }x/{\nu }_f$ is the local Reynolds number. Note that, the reduced Sherwood number follows the reduced Nusselt number due to the passive control condition in Equation (13) for $Nb=Nt$. The effect of the changes in the boundary condition is instantly obvious where the parameters $Nb$ and $Nt$ now appear in the boundary conditions as well as in the differential equations.

3. Stability Analysis

The existence of the non-uniqueness solutions of Equations (10) to (13) is observed for a certain range of the physical parameters. A temporal stability analysis is therefore needed to determine which solution is stable and thus physically reliable in the long run. This technique was initiated by Merkin [38] in 1986. A dimensionless time variable $\tau$ was introduced by Weidman et al. [39] to further study the stability of the solutions as time passes. They concluded that the first solutions are stable, while the second solutions are unstable. As in Weidman et al. [39], the new variables based on Equation (6) are given as

$$\psi ={\nu }_{f}xf\left(\eta ,\tau \right), \theta \left(\eta ,\tau \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }}, \varphi \left(\eta ,\tau \right)=\frac{C-C_{\infty }}{C_{\infty }}, \eta =\frac{U_{\infty } r^2}{{\nu }_f x}, \tau =\frac{2U_{\infty }t}{x}$$

(18)

The unsteady form of Equations (1) to (4) are considered to analyze the stability of the solutions. Using (18) and following the similar approach as in Section 2, we get

$$2\frac{{\mu }_{hnf}/{\mu }_f}{{\rho }_{hnf}/{\rho }_f}\frac{\partial }{\partial \eta }\left(\eta \frac{{\partial }^{2}f}{\partial {\eta }^2}\right)+f\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{{\partial }^{2}f}{\partial \eta \partial \tau }-\tau \left(\frac{\partial f}{\partial \tau }\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{\partial f}{\partial \eta }\frac{{\partial }^{2}f}{\partial \eta \partial \tau }\right)=0$$

(19)

$$\frac{2}{\mathrm{Pr}}\frac{k_{hnf}/k_f}{{(\rho C_p)}_{hnf}/{(\rho C_p)}_f}\frac{\partial }{\partial \eta }\left(\eta \frac{\partial \theta }{\partial \eta }\right)+f\frac{\partial \theta }{\partial \eta }+2\eta \left(Nb\frac{\partial \varphi }{\partial \eta }\frac{\partial \theta }{\partial \eta }+Nt{\left(\frac{\partial \theta }{\partial \eta }\right)}^2\right)-\frac{\partial \theta }{\partial \tau }-\tau \left(\frac{\partial f}{\partial \tau }\frac{\partial \theta }{\partial \eta }-\frac{\partial f}{\partial \eta }\frac{\partial \theta }{\partial \tau }\right)=0$$

(20)

$$2\frac{\partial }{\partial \eta }\left(\eta \frac{\partial \varphi }{\partial \eta }\right)+Scf\frac{\partial \varphi }{\partial \eta }+2\frac{Nt}{Nb}\frac{\partial }{\partial \eta }\left(\eta \frac{\partial \theta }{\partial \eta }\right)-Sc\frac{\partial \varphi }{\partial \tau }-Sc\tau \left(\frac{\partial f}{\partial \tau }\frac{\partial \varphi }{\partial \eta }-\frac{\partial f}{\partial \eta }\frac{\partial \varphi }{\partial \tau }\right)=0$$

(21)

subject to:

$$\begin{array}{c}f\left(c,\tau \right)-\tau \frac{\partial f}{\partial \tau }\left(c,\tau \right)=\lambda c, \frac{\partial f}{\partial \eta }\left(c,\tau \right)=\frac{\lambda }{2}, \theta \left(c,\tau \right)=1, Nb\frac{\partial \varphi }{\partial \eta }\left(c,\tau \right)+Nt\frac{\partial \theta }{\partial \eta }\left(c,\tau \right)=0\\ \frac{\partial f}{\partial \eta }\left(\eta ,\tau \right)\to \frac{1}{2}, \theta \left(\eta ,\tau \right)\to 0, \varphi \left(\eta ,\tau \right)\to 0 \mathrm{as} \eta \to \infty .\end{array}$$

(22)

To examine the stability behavior, the disturbance is imposed to the steady solution $f\left(\eta \right)=f_0\left(\eta \right)$, $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$ and $\varphi \left(\eta \right)={\varphi }_0\left(\eta \right)$ of Equations (10) to (13) by using the following relations (see [39]):

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

(23)

where $\gamma$ indicates the unknown eigenvalue that determines the stability of the solutions and $F\left(\eta \right)$, $G\left(\eta \right)$ and $H\left(\eta \right)$ are comparatively small to $f_0\left(\eta \right)$, ${\theta }_0\left(\eta \right)$ and ${\varphi }_0\left(\eta \right)$. The disturbance is taken exponentially as it demonstrates the rapid decline or development of the disturbance. By inserting Equation (23) into Equations (19) to (21) and by setting $\tau =0$, we obtain

$$2\frac{{\mu }_{hnf}/{\mu }_f}{{\rho }_{hnf}/{\rho }_f}{\left(\eta F^{″}\right)}^{\prime }+f_0{F}^{″}+f_0^{″}F+\gamma F^{\prime }=0$$

(24)

$$\begin{array}{c}\frac{2}{\mathrm{Pr}}\frac{k_{hnf}/k_f}{{(\rho C_p)}_{hnf}/{(\rho C_p)}_f}{\left(\eta G^{\prime }\right)}^{\prime }+f_0{G}^{\prime }+{\theta }_0^{\prime }F\\ +2\eta \left(Nb{\varphi }_0^{\prime }{G}^{\prime }+Nb{\theta }_0^{\prime }{H}^{\prime }+2Nt{\theta }_0^{\prime }{G}^{\prime }\right)+\gamma G=0\end{array}$$

(25)

$$2{\left(\eta H^{\prime }\right)}^{\prime }+Sc{f}_0{H}^{\prime }+Sc{\varphi }_0^{\prime }F+2\frac{Nt}{Nb}{\left(\eta G^{\prime }\right)}^{\prime }+\gamma ScH=0$$

(26)

subject to:

$$\begin{array}{c}F\left(c\right)=0, F^{\prime }\left(c\right)=0, G\left(c\right)=0, Nb{H}^{\prime }\left(c\right)+Nt{G}^{\prime }\left(c\right)=0\\ F^{\prime }\left(\eta \right)\to 0, G\left(\eta \right)\to 0, H\left(\eta \right)\to 0 \mathrm{as} \eta \to \infty \end{array}$$

(27)

Without loss of generality, the values of $\gamma$ from Equations (24) to (27) are obtained for the case of $F^{″}\left(c\right)=$ 1 as suggested and discussed by Harris et al. [40].

4. Numerical Method

The bvp4c solver in Matlab software (Matlab_R2014b, MathWorks, Singapore) is utilized for evaluating the Equations (10) to (13), numerically. As described in Shampine et al. [41], the aforesaid solver occupies a finite difference method that employs the 3-stage Lobatto IIIa formula. The selection of the initial guess and the boundary layer thickness, ${\eta }_{\infty }$ are particular depend on the parameters applied to obtain the solutions. Moreover, several researchers [42,43,44,45,46] are also employing this solver for solving the boundary layer flow problems. First, Equations (10) to (12) are reduced to a system of ordinary differential equations of the first order. Now, Equation (10) may be written as

$$\begin{array}{r}f=y\left(1\right)\\ f^{\prime }=y^{\prime }\left(1\right)=y\left(2\right)\end{array}$$

(28a)

$$f^{″}=y^{\prime }\left(2\right)=y\left(3\right)$$

(28b)

$$f^{‴}=y^{\prime }\left(3\right)=-\frac{1}{2\eta }\left\{\frac{{\rho }_{hnf}/{\rho }_f}{{\mu }_{hnf}/{\mu }_f}y\left(1\right)y\left(3\right)+2y\left(3\right)\right\}$$

(28c)

while Equation (11) reduces to:

$$\begin{array}{r}\theta =y\left(4\right)\\ {\theta }^{\prime }=y^{\prime }\left(4\right)=y\left(5\right)\end{array}$$

(29a)

$${\theta }^{″}=y^{\prime }\left(5\right)=-\frac{1}{2\eta }\left\{\mathrm{Pr}\frac{{\left(\rho C_p\right)}_{hnf}/{\left(\rho C_p\right)}_f}{k_{hnf}/k_f}\left[y\left(1\right)y\left(5\right)+2\eta \left(Nby\left(7\right)y\left(5\right)+Nty{\left(5\right)}^2\right)\right]+2y\left(5\right)\right\}$$

(29b)

and Equation (12) reduces to:

$$\begin{array}{r}\varphi =y\left(6\right)\\ {\varphi }^{\prime }=y^{\prime }\left(6\right)=y\left(7\right)\end{array}$$

(30a)

$${\varphi }^{″}=y^{\prime }\left(7\right)=-\frac{1}{2\eta }\left\{Scy\left(1\right)y\left(7\right)+2\frac{Nt}{Nb}{\left(\eta y\left(5\right)\right)}^{\prime }+2y\left(7\right)\right\}$$

(30b)

with the boundary conditions:

$$\begin{array}{c}ya\left(1\right)=\lambda c, ya\left(2\right)=\frac{\lambda }{2}, ya\left(4\right)=1, Nbya\left(7\right)+Ntya\left(5\right)=0\\ yb\left(2\right)\to \frac{1}{2}, yb\left(4\right)\to 0, yb\left(6\right)\to 0\end{array}$$

(31)

Then, Equations (28) to (31) are coded in Matlab software (Matlab_R2014b, MathWorks, Singapore) and their solutions are obtained by the bvp4c solver. The solver will then run, and the outcomes will be printed out as numerical solutions and graphs.

5. Results and Discussion

In the present study, we consider various volume fractions of Cu $\left({\phi }_2\right)$, while the volume fraction of Al₂O₃ is kept fixed at ${\phi }_1=0.1$ and water as the base fluid. To ensure the accuracy of the computation, the present results are validated with the existing data from the previous studies. The values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ when ${\phi }_1={\phi }_2=0$ (regular fluid), $\lambda =\mathrm{Sc}=0$ and $\mathrm{Pr}=0.733$ for various values of $c$ in the absence of $Nt$ and $Nb$ are displayed in

Table 3. Values of $f^{″}\left(c\right)$ for regular fluid (${\phi }_1={\phi }_2=0$) when $\lambda =0$ with various values of $c$.

$\mathit{c}$	Chen and Smith [26]	Ishak et al. [30]	Grosan and Pop [32]	Soid et al. [33]	Present Results
0.1	1.28881	1.2888	1.289074	1.288778	1.288778
0.01	8.49244	8.4924	8.492173	8.491454	8.491454
0.001	62.16372	62.1637	62.161171		62.158227

and

Table 4. Values of $-{\theta }^{\prime }\left(c\right)$ for regular fluid (${\phi }_1={\phi }_2=0$) when $\lambda =\mathrm{Sc}=0$ and $\mathrm{Pr}=0.733$ for various values of $c$ in the absence of $Nt$ and $Nb$.

$\mathit{c}$	Chen and Smith [26]	Grosan and Pop [32]	Present Results
0.1	2.434	2.441675	2.439692
0.01		16.306544	16.283107
0.001		120.55034	120.264815

, respectively. The results are comparable with those obtained by the previous studies. On the other hand,

Table 5. Values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ for Cu/water nanofluid when ${\phi }_1=\lambda =\mathrm{Sc}=0$ and $\mathrm{Pr}=7$ for various values of $c$ and ${\phi }_2$ in the absence of $Nt$ and $Nb$.

		${\mathit{f}}^{″}\left(\mathit{c}\right)$		$-{\mathit{\theta }}^{\prime }\left(\mathit{c}\right)$
$\mathit{c}$	${\mathit{\phi }}_2$	Grosan and Pop [32]	Present Results	Grosan and Pop [32]	Present Results
0.1	0.05	1.347208	1.347125	3.682009	3.681817
	0.1	1.382008	1.381635	3.586544	3.586427
	0.2	1.404136	1.404050	3.389762	3.389682
0.01	0.05	8.771680	8.771503	22.284916	22.284751
	0.1	8.935933	8.935143	21.816182	21.816075
	0.2	9.041011	9.040694	20.877668	20.877591
0.001	0.05	63.884384	63.718195	153.570113	153.569704
	0.1	64.653616	64.621235	150.977665	150.977406
	0.2	65.235057	65.200519	145.860889	145.860701

is provided to describe the values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ for Cu/water nanofluid when ${\phi }_1=\lambda =\mathrm{Sc}=0$ and $\mathrm{Pr}=7$ with various values of $c$ and ${\phi }_2$ in the absence of $Nt$ and $Nb$. The present numerical results are compared with the previous results obtained by Grosan and Pop [32] and show a favorable agreement. The increment values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ are observed for the smaller needle size as can be seen in Table 3, Table 4 and Table 5. It is also observed that the increasing values of ${\phi }_2$ tend to increase the values of $f^{″}\left(c\right)$ but decrease the values of $-{\theta }^{\prime }\left(c\right)$ as shown in Table 5.

The plots of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ for selected parameters with various values of ${\phi }_2$ and $c$ against $\lambda$ are portrayed in Figure 2, Figure 3, Figure 4 and Figure 5, respectively. Results show that dual solutions exist for a certain range of $\lambda$. We observe that the critical values ${\lambda }_c$ increase and move slightly to the right with the increasing of ${\phi }_2$ and $c$. From our computation, ${\lambda }_c=-1.13048$, $-1.10314$ and $-1.08436$ are the critical values for ${\phi }_2=0.01$, $0.05$ and $0.1$, respectively as shown in Figure 2 and Figure 3. Furthermore, the increasing of ${\phi }_2$ has the tendency to decrease the skin friction coefficient for $\lambda <-0.7$ but increase the heat transfer rate for $\lambda <-0.3$. We also observe that the heat transfer rate is almost the same for $\lambda =-0.3$. Meanwhile, the critical values of $\lambda$ for $c=0.1$ and $0.2$ are ${\lambda }_c=-1.08436$ and $-0.90271$ as can be seen in Figure 4 and Figure 5. It is observed that the values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ are greater for smaller $c$.

The impact of several parameters on the velocity $f^{\prime }\left(\eta \right)$, temperature $\theta \left(\eta \right)$ and nanoparticle concentration $\varphi \left(\eta \right)$ profiles are presented in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. These profiles asymptotically satisfy the free stream conditions (13), thus giving us confidence in the accuracy of the numerical solutions. Figure 6, Figure 7 and Figure 8 depict the effect of ${\phi }_2$ on velocity profiles $f^{\prime }\left(\eta \right)$, temperature profiles $\theta \left(\eta \right)$ and concentration profiles $\varphi \left(\eta \right)$ when $\lambda =-1$, $Sc=Nb=Nt=0.5$, $c={\phi }_1=0.1$ and $\mathrm{Pr}=6.2$. Small effects are observed for the first solutions (upper branch) of $f^{\prime }\left(\eta \right)$, $\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ with rising ${\phi }_2$, but the effects are more pronounced for the second solutions (lower branch).

The profiles of $f^{\prime }\left(\eta \right)$, $\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ for different values of $\lambda$ when $Sc=Nb=Nt=0.5$, $c={\phi }_1={\phi }_2=0.1$ and $\mathrm{Pr}=6.2$ are plotted in Figure 9, Figure 10 and Figure 11, respectively. It is observed that the first solution approaching the second solution as the critical value ${\lambda }_c$ is approached. This observation is consistent with the results presented in Figure 2, Figure 3, Figure 4 and Figure 5. Moreover, the effect of $Nt$ on $\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ when $\lambda =-1$, $Sc=Nb=0.5$, $c={\phi }_1={\phi }_2=0.1$ and $\mathrm{Pr}=6.2$ are displayed in Figure 12 and Figure 13. It is clear that $\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ decrease for both branches as $Nt$ increases. Figure 14 exhibits the behavior of $\varphi \left(\eta \right)$ with the variation of $Nb$ when $\lambda =-1$, $Sc=Nt=0.5$, $c={\phi }_1={\phi }_2=0.1$ and $\mathrm{Pr}=6.2$. A rise in $Nb$ produces an increase of $\varphi \left(\eta \right)$ for both branches. Obviously, $\theta \left(\eta \right)$ does not affected by $Nb$ due to the passive control condition applied in (13) and consequently, the energy Equation (11) does not contain the parameter $Nb$.

Figure 15 displays the smallest eigenvalue $\gamma$ against $\lambda$ when $c={\phi }_1={\phi }_2=0.1$. As described in Equation (23), the flow is stable when there is an initial decay of disturbance as time passes, i.e., $e^{-\gamma \tau }\to 0$ as $\tau \to \infty$. This will happen for $\gamma >0$. Meanwhile, the flow is unstable for $\gamma <0$ due to the initial growth of disturbance as $\tau \to \infty$. From Figure 15, we notice that the values of $\gamma$ are positive for the first solutions (upper branch), while they are negative for the second solutions (lower branch). Also, the values of $\gamma$ approach to zero for both branches when $\lambda \to {\lambda }_c=-1.08436$. Thus, this finding confirms that the first solutions are stable and physically reliable in the long run, while the second solutions are not. Furthermore, we also conclude that the bifurcation of the solutions happens at the critical value $\lambda ={\lambda }_c$.

6. Conclusions

In this paper, the flow past a permeable moving thin needle in a hybrid nanofluid with thermophoresis and Brownian motion effects was studied. The results validation was done for the limiting cases where the present results are comparable with the existing results. The present results revealed that the added hybrid nanoparticles led to the increment of the heat transfer rate. Dual solutions were obtained for a certain range of the constant moving parameter $\lambda$. Furthermore, it was discovered that the bifurcation of the solutions occurs in the region $\lambda <0$, i.e., when the needle moved toward the origin. The critical values of $\lambda$ slightly increased for larger values of ${\phi }_2$ and $c$. Furthermore, the values of $f^{″}\left(c\right)$ and $-{\theta }^{\prime }\left(c\right)$ are higher for smaller $c$. It was also observed that $\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ decreased for larger values of $Nt$. The increase of $Nb$ led to an increase in $\varphi \left(\eta \right)$. It was shown that between the two solutions, only one of them is stable, while the other is unstable as time evolved.