Stability Analysis of Stagnation-Point Flow in a Nanofluid over a Stretching/Shrinking Sheet with Second-Order Slip, Soret and Dufour Effects: A Revised Model

Abstract

The mathematical model of the two-dimensional steady stagnation-point flow over a stretching or shrinking sheet of nanofluid in the presence of the Soret and Dufour effects and of second-order slip at the boundary was considered in this paper. The partial differential equations were transformed into the ordinary differential equations by applying a suitable similarity transformation. The numerical results were obtained by using bvp4c codes in Matlab. The skin friction coefficient, heat transfer coefficient, mass transfer coefficient, as well as the velocity, temperature, and concentration profiles were presented graphically for different values of slip parameters, Soret effect, Dufour effect, Brownian motion parameter, and thermophoresis parameter. A dual solution was obtained in this present paper. The presence of the slip parameters (first- and second-order slip parameters) was found to expand the range of solutions. However, the presence of the slip parameters led to a decrease in the skin friction coefficient, whereas the heat transfer coefficient increased. Besides that, a larger Soret effect (smallest Dufour effect) led to the decrement of the heat transfer coefficient. The effects of the Brownian motion and thermophoresis parameters on the heat transfer coefficient were also studied in this paper. A stability analysis was performed in this paper to verify the stability of the solutions obtained.

1. Introduction

A nanofluid can be defined as the dispersion of nanoparticles into a base fluid in which the collision between the nanoparticles would enhance the thermal conductivity of the fluid (see Masuda et al. [1]). By adding or suspending nanoparticles into a basic fluid (water, oil, and ethylene glycol), the effectiveness of the thermal conductivity of the nanofluid is expected to increase twice because of the larger nanoparticle surface area. Thus, nanofluids have a great potential in heat transfer applications. The applications of nanofluids are numerous in the industrial and automotive sector as well as in electronic devices. For instance, the use of nanofluids can be found in refrigeration, lubrication, radiators, heat exchangers, heat transfer, and also in the cooling of microchips devices in laptops and cell phones. The two models that have been introduced to study the characteristics of nanofluids are the Tiwari and Das model [2] and the Buongiorno model [3]. The main objective of the Tiwari and Das model is to study the effects of nanoparticles volumetric fraction, as reported by Ul Haq et al. [4]. On the other hand, the Buongiorno model focuses on the effects of Brownian motion and thermophoresis in nanofluids. The proposed model has been implemented by Nield and Kuznetsov [5], Kuznetsov and Nield [6], Bachok et al. [7], Mustafa et al. [8], Hamid et al. [9], Mansur et al. [10], and Jamaludin et al. [11]. Later, the boundary condition of the Buongiorno model was modified by Kuznetsov and Nield [12,13] resulting in a revised model. These researchers improved the Buongiorno model by extending it to the case where the nanofluid particles fraction on the boundary is controlled passively rather than actively. Zaimi et al. [14] reported the flow past a stretching/shrinking sheet with suction in a nanofluid using the revised model. Different approaches have been made by Esfe et al. [15] to present a numerical simulation of the natural convection in an enclosure filled with nanofluid when the cylindrical block is heated. Furthermore, several works on solving nanofluid problems using the lattice Boltzmann method were reported by Karimipour [16] and Karimipour et al. [17,18]. Esfe et al. [19] experimentally investigated the effects of the solid volume fraction and thermal conductivity of nanofluids.

The consideration of the slip conditions at the boundary in certain cases should not be ignored, such as in rarefied gas flows in micro-scale devices and low-pressure situations [20,21,22]. The new formulation of the second-order slip flow was provided by Wu [23] using the kinetic theory. The new slip model has then been used by Fang et al. [24] and Fang and Aziz [25]. Nandeppanavar et al. [26] analyzed the Wu’s model over a stretching sheet with no-linear Navier boundary condition. The second-order slip flow and heat transfer over a permeable shrinking surface embedded in a porous medium and in mixed convection flow were reported by Yasin et al. [27] and Singh and Chamkha [28], respectively. The work by Khader et al. [29] described the second-order slip embedded in a porous medium when the viscous dissipation was taken into account over a permeable stretching sheet. Hakeem et al. [30] studied the effect of a magnetic field on second- order slip flow in a nanofluid over a stretching/shrinking sheet with radiation. Soid et al. [31] applied the new slip model and temperature jump to an axisymmetric stagnation-point flow over a stretching/shrinking surface. The application in industrial processes involving slip velocity can be found in polymer solutions and emulsions process in the oil industry.

According to Platten [32], the Swiss scientist named Charles Soret introduced in 1879 the Soret effect, also known as thermodiffusion, which can be defined as a diffusion of molecules or particles from a higher temperature towards a lower temperature resulting from mass flux. Besides that, the Dufour effect or diffusion-thermo can be defined as a diffusion of molecules or particles from a higher concentration to a lower concentration as a result of energy flux. The application of the Soret and Dufour effects can be found in the fabrication of semiconductor devices in molten metal, semiconductor mixtures, separation of polymers and DNA, as well as in optimum oil recovery from hydrocarbon reservoirs [33]. An extensive literature discussing the Soret and Dufour effects exists [34,35,36,37,38,39,40,41].

A few years back, the stability analysis gained attention among researchers because of the presence of dual or more solutions in their work. The purpose of performing the stability analysis was to determine or verify if a solution was stable and physically realizable. The first pioneer who proposed the stability solution in his study was Merkin [42]. The implementation of the stability analysis was done by Merill et al. [43], Weidman et al. [44], and Harris et al. [45] in their works. After that, Mahapatra and Nandy [46] considered the stability analysis over a power law shrinking surface. The consideration of stability solutions in nanofluids over a stretching/shrinking sheet was also reported by Noor et al. [47], Nazar et al. [48], and Jusoh and Nazar [49]. The stability analysis of an unsteady three-dimensional boundary layer flow with suction was solved by Hafidzuddin et al. [50]. Ishak [51] and Ismail et al. [52] performed a stability analysis over a linear and exponentially shrinking surface, respectively. A stability analysis over a moving plate was reported by Najib et al. [53,54], Bachok et al. [55], and also Noor et al. [56].

Motivated by the above studies, the present paper aims to investigate the effects of Soret and Dufour by following the governing equations proposed by Bhattacharyya et al. [36] on the stagnation-point flow of a nanofluid over a stretching/shrinking surface with second-order slip conditions (following the formulated second-order slip model by Wu [23]) at the boundary with the improvised model proposed by Kuznetsov and Nield [12]. The governing equations in the form of partial differential equations were transformed into ordinary differential equations using suitable similarity variables and then solved numerically using a solver in Matlab (Matlab R2013a, Mathwork, Natick, MA, USA, 1984). The stability solutions were analyzed to determine the stability of the solutions obtained.

2. Basic Equations

We considered the flow of an incompressible nanofluid in the region $y>0$ driven by a stretching/shrinking surface located at $y=0$ with a fixed stagnation point at $x=0$ in the presence of Soret and Dufour effects as well as of a second-order slip effect at the boundary, as illustrated in Figure 1. It was assumed that the free stream and stretching/shrinking velocity $\frac{U_w}{U_{\infty }}$ were in the linear form, with $U_{\infty }=cx$ and $U_w=bx$, respectively, where $c$ and $b$ were constant with $c>0$. It should be noted that $b>0$ and $b<0$ correspond to the stretching and shrinking sheet, respectively. Under these conditions, the boundary layer equations are

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

(1)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=U_{\infty }\frac{d{U}_{\infty }}{dx}+\nu \frac{{\partial }^{2}u}{\partial y^2}$$

(2)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2}+\mathsf{\Omega }\left(D_B\frac{\partial C}{\partial y}\frac{\partial T}{\partial y}+\frac{D_T}{T_{\infty }}{\left(\frac{\partial T}{\partial y}\right)}^2\right)+\frac{D_B{k}_T}{c_s{c}_p}\frac{{\partial }^{2}C}{\partial y^2}$$

(3)

$$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_B\frac{{\partial }^{2}C}{\partial y^2}+\frac{D_T}{T_{\infty }}\frac{{\partial }^{2}T}{\partial y^2}+\frac{D_B{k}_T}{T_m}\frac{{\partial }^{2}T}{\partial y^2}$$

(4)

along with the initial and boundary conditions

$$\begin{array}{ccc}\hfill t<0:& u=v=0,\hspace{1em}T=T_{\infty },\hspace{1em}C=C_{\infty }\hfill & \mathrm{for}\mathrm{any}x,y\hfill \\ \hfill t\ge 0:& u=U_w+A{\left(\frac{\partial u}{\partial y}\right)}_{y=0}+B{\left(\frac{{\partial }^{2}u}{\partial y^2}\right)}_{y=0},\hspace{1em}v=0,\hfill & \\ & T=T_w,\hspace{1em}{D}_B\frac{\partial C}{\partial y}+\frac{D_T}{T_{\infty }}\frac{\partial T}{\partial y}=0\hfill & \mathrm{at}y=0\hfill \\ & u\to U_{\infty },\hspace{1em}T\to T_{\infty },\hspace{1em}C\to C_{\infty }\hfill & \mathrm{as}y\to \infty \hfill \end{array}$$

(5)

where $u$ and $v$ are the velocity components along $x-$ and $y-$ axes, respectively, $T$ is the nanofluid temperature. $U_{slip}$ is the slip velocity at the wall. The Wu’s slip velocity model (valid for arbitrary Knudsen numbers, Kn) is used in this paper and is given as follows [23]

$$U_{slip}=A\frac{\partial u}{\partial y}+B\frac{{\partial }^{2}u}{\partial y^2}$$

(6)

where $A$ and $B$ are constant.

3. Steady-State Solution $\mathbf{\left(}\mathbf{\partial }\mathbf{/}\mathbf{\partial }\mathbf{t}\mathbf{=}\mathbf{0}\mathbf{\right)}$

The following similarity transformation was introduced

$$\psi =\sqrt{c\nu }xf(\eta ),\hspace{1em}\theta (\eta )=\frac{T-T_{\infty }}{T_w-T_{\infty }},\hspace{1em}\varphi (\eta )=\frac{C-C_{\infty }}{C_w-C_{\infty }},\hspace{1em}\eta =\sqrt{\frac{c}{\nu }}y,$$

(7)

where $\eta$ is the similarity variable, $\psi$ is the stream function defined as $u=\partial \psi /\partial y$, and $v=- \partial \psi /\partial x$, which automatically satisfies Equation (1). The similarity variables (7) were substituted by Equations (2)–(4) to obtain the following ordinary (similarity) differential equations

$$f^{‴}+f{f}^{″}-{f^{\prime }}^2+1=0$$

(8)

$$\frac{1}{\mathrm{Pr}}{\theta }^{″}+f{\theta }^{\prime }+Nb{\theta }^{\prime }{\varphi }^{\prime }+Nt{{\theta }^{\prime }}^2+Du{\varphi }^{″}=0$$

(9)

$${\varphi }^{″}+Lef{\varphi }^{\prime }+\frac{Nt}{Nb}{\theta }^{″}+LeSr{\theta }^{″}=0$$

(10)

subject to the boundary conditions

$$\begin{array}{llll}f\left(0\right)=0,& f^{\prime }\left(0\right)=\epsilon +\sigma f^{″}\left(0\right)+\delta f^{‴}\left(0\right),& \theta \left(0\right)=1,& Nb{\varphi }^{\prime }\left(0\right)+Nt{\theta }^{\prime }\left(0\right)=0,\\ f^{\prime }(\infty )\to 1,& \theta (\infty )\to 0,\hspace{1em}\varphi (\infty )\to 0& & \end{array}$$

(11)

In the above equations, primes denote the differentiation with respect to $\eta$. According to Mukhopadhyay and Andersson [57], $A=\sigma \sqrt{\nu /c}$ and $B=\delta (\nu /c)$, with $\sigma >0$ being the first velocity slip parameter, and $\delta <0$ the second velocity slip parameter (see Fang et al. [24]). Here, $\mathrm{Pr}$ is the Prandt number, $Le$ is the Lewis number, $Sr$ is the Soret number, $Du$ is the Dufour number, $Nb$ is the Brownian motion, $Nt$ is the thermophoresis parameter, and $\epsilon$ is the velocity ratio parameter, which are defined as:

$$\begin{array}{l}\mathrm{Pr}=\frac{\nu }{\alpha },\hspace{1em}\mathrm{Le}=\frac{\nu }{D_B},\hspace{1em}Sr=\frac{D_B{K}_T}{\nu T_m}\frac{T_w-T_{\infty }}{C_w-C_{\infty }},\hspace{1em}Du=\frac{D_B{K}_T}{\nu c_s{c}_p}\frac{C_w-C_{\infty }}{T_w-T_{\infty }},\\ \epsilon =\frac{b}{c},\hspace{1em}Nb=\frac{\mathsf{\Omega }{D}_B}{\nu }\left(C_w-C_{\infty }\right),\hspace{1em}Nt=\frac{\mathsf{\Omega }{D}_T}{T_{\infty }\nu }\left(T_w-T_{\infty }\right)\end{array}$$

(12)

where $\epsilon >0$ for stretching and $\epsilon <0$ for shrinking.

The physical quantities of practical interest are the local 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 U_{\infty }{}^2},\hspace{1em}N{u}_x=\frac{x{q}_w}{k(T_w-T_{\infty })},\hspace{1em}S{h}_x=\frac{x{q}_m}{D_B(C_w-C_{\infty })},$$

(13)

where ${\tau }_w$ is the skin friction or the shear stress on the stretching/shrinking sheet, $q_w$ is the heat flux from the surface of the plate, and $q_m$ is the mass flux from the surface of the plate, which are given by

$${\tau }_w=\mu {\left(\frac{\partial u}{\partial y}\right)}_{y=0},\hspace{1em}{q}_w=-k{\left(\frac{\partial T}{\partial y}\right)}_{y=0},\hspace{1em}{q}_m=-D_B{\left(\frac{\partial C}{\partial y}\right)}_{y=0},$$

(14)

Using (7) in (13) and (14), we obtained

$${\mathrm{Re}}_x{}^{1/2}{C}_f=f^{″}(0),\hspace{1em}{\mathrm{Re}}_x{}^{-1/2}N{u}_x=-{\theta }^{\prime }(0),\hspace{1em}{\mathrm{Re}}_x{}^{-1/2}S{h}_x=-{\varphi }^{\prime }(0)$$

(15)

where ${\mathrm{Re}}_x=c{x}^2/\nu$ is the local Reynolds number.

4. Stability Analysis

Weidman et al. [44] and Roşca and Pop [58] have shown that the lower branch solutions (second solutions) are unstable (not realizable physically), while the upper branch solutions (first solutions) are stable (physically realizable) by considering the unsteady Equations (2)–(4). Thus, the new dimensionless time variable $\tau =ct$ was introduced. The use of $\tau$ was associated with an initial value problem and was consistent with the question of which solution will be obtained in practice (physically realizable). Using the variables $\tau$ and (7), we have

$$\begin{array}{l}\psi =\sqrt{c\nu }xf(\eta ,\tau ),\hspace{1em}\theta (\eta ,\tau )=\frac{T-T_{\infty }}{T_w-T_{\infty }},\hspace{1em}\varphi (\eta ,\tau )=\frac{C-C_{\infty }}{C_w-C_{\infty }},\\ \eta =\sqrt{\frac{c}{\nu }}y,\hspace{1em}\tau =ct\end{array}$$

(16)

Equations (2) and (3) can be written as

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

(17)

$$\frac{1}{\mathrm{Pr}}\frac{{\partial }^2\theta }{\partial {\eta }^2}+f\frac{\partial \theta }{\partial \eta }+Nb\frac{\partial \theta }{\partial \eta }\frac{\partial \varphi }{\partial \eta }+Nt{\left(\frac{\partial \theta }{\partial \eta }\right)}^2+Du\frac{{\partial }^2\varphi }{\partial {\eta }^2}-\frac{\partial \theta }{\partial \tau }=0$$

(18)

$$\frac{{\partial }^2\varphi }{\partial {\eta }^2}+Lef\frac{\partial \varphi }{\partial \eta }+\frac{Nt}{Nb}\frac{{\partial }^2\theta }{\partial {\eta }^2}+LeSr\frac{{\partial }^2\theta }{\partial {\eta }^2}-Le\frac{\partial \varphi }{\partial \tau }=0$$

(19)

subject to the boundary conditions

$$\begin{array}{l}f(0,\tau )=0,\hspace{1em}\frac{\partial f}{\partial \eta }(0,\tau )=\epsilon +\sigma \frac{{\partial }^{2}f}{\partial {\eta }^2}+\delta \frac{{\partial }^{3}f}{\partial {\eta }^3},\\ \theta (0,\tau )=1,\hspace{1em}Nb\frac{\partial \varphi }{\partial \eta }(0,\tau )+Nt\frac{\partial \theta }{\partial \eta }(0,\tau )=0,\\ \frac{\partial f}{\partial \eta }(\eta ,\tau )\to 1,\hspace{1em}\theta (\eta ,\tau )\to 0,\hspace{1em}\varphi (\eta ,\tau )=0\hspace{1em}\mathrm{as}\hspace{1em}\eta \to \infty \end{array}$$

(20)

To determine the stability of the solution, with $f=f_0(\eta )$, $\theta ={\theta }_0(\eta )$, and $\varphi ={\varphi }_0(\eta )$ satisfying the boundary value problem (17)–(20), we have (see Weidman et al. [44] or Roşca and Pop, [58])

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

(21)

where $\gamma$ is an unknown eigenvalue parameter and $F(\eta )$, $G(\eta )$ and $H(\eta )$ are small relative to $f=f_0(\eta )$, $\theta ={\theta }_0(\eta )$, and $\varphi ={\varphi }_0(\eta )$. The solution of the eigenvalue problem (17)–(20) gives an infinite set of eigenvalues ${\gamma }_1<{\gamma }_2<{\gamma }_3\dots$; if ${\gamma }_1$ is negative, there is an initial growth of disturbances, and the flow is unstable, but when ${\gamma }_1$ is positive, there is an initial decay, and the flow is stable. Introducing (21) into (17)–(20), we obtain the following linearized problem

$${F_0}^{‴}+f_0{F_0}^{″}+{f_0}^{″}{F}_0-2{f_0}^{\prime }{F_0}^{\prime }+\gamma {F_0}^{\prime }=0$$

(22)

$$\frac{1}{\mathrm{Pr}}{G_0}^{″}+f_0{G_0}^{\prime }+F_0{{\theta }_0}^{\prime }+Nb{{\theta }_0}^{\prime }{H_0}^{\prime }+Nb{G_0}^{\prime }{{\varphi }_0}^{\prime }+2Nt{{\theta }_0}^{\prime }{G_0}^{\prime }+Du{H_0}^{″}+\gamma G_0=0$$

(23)

$${H_0}^{″}+Le{f}_0{H_0}^{\prime }+Le{F}_0{{\varphi }_0}^{\prime }+\frac{Nt}{Nb}{G_0}^{″}+LeSr{G_0}^{″}+Le\gamma H_0=0$$

(24)

subject to the boundary conditions

$$\begin{array}{l}{F}_0(0)=0,\hspace{1em}{F_0}^{\prime }(0)=\sigma {F_0}^{″}(0)+\delta {F_0}^{‴}(0),\\ G_0(0)=0,\hspace{1em}Nb{H_0}^{\prime }\left(0\right)+Nt{G_0}^{\prime }\left(0\right)=0,\\ {F_0}^{\prime }(\eta )\to 0,\hspace{1em}{G}_0(\eta )\to 0,\hspace{1em}{H}_0(\eta )\to 0\hspace{1em}\mathrm{as}\hspace{1em}\eta \to \infty \end{array}$$

(25)

It is inevitable to mention that for particular values of $\epsilon , \sigma , \delta , \mathrm{Pr}, Sr, Du, Nb, Nt$, and $Le$, the stability of the corresponding steady flow solution $f_0(\eta )$, ${\theta }_0(\eta )$, and ${\varphi }_0(\eta )$, was determined by the smallest eigenvalue $\gamma$. According to Harris et al. [45], the range of possible eigenvalues can be determined by relaxing a boundary condition either on $F_0(\eta )$, $G_0(\eta )$ or $H_0(\eta )$. For the present problem, we relaxed the condition that ${F_0}^{\prime }(\eta )\to 0$, as $\eta \to \infty$ and, for a fixed value of $\gamma$, we solved the system (22)–(24) subject to (25) considering the new boundary condition ${F_0}^{″}(0)=1$.

5. Results and Discussion

To discuss in detail the effects of the slip parameters (first-order slip $\sigma$ and second-order slip $\delta$), the Soret and Dufour effects ($Sr$ and $Du$), the Brownian motion $Nb$, and also the thermophoresis parameter $Nt$, a set of ordinary differential Equations (8)–(10) along with boundary conditions (11) was substituted into bvp4c codes in Matlab software to obtain the numerical results of the tested parameters. The results of the skin friction coefficient, heat and mass transfer coefficient, as well as the velocity, temperature, and concentration profiles were graphically presented (Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10). The effects of the first-order slip parameter $\sigma$ and second-order slip parameter $\delta$ were plotted (Figure 2, Figure 3 and Figure 4). From the figures, it was clear that a unique solution existed for $\epsilon >-1$, dual solutions were found for ${\epsilon }_c\le \epsilon \le -1$, and no solution was reported for $\epsilon <{\epsilon }_c$. These three figures indicate that an increase in the slip parameters ($\sigma$ and $\delta$) caused a decrement in the value of the skin friction coefficient. On the other hand, an increment of the heat transfer coefficient was observed for an increase of the slip parameters. A different result can be observed in Figure 3 and Figure 4, in which the range of solutions for Figure 3 was widely expanded for the increment of the first-order slip parameter, while, in Figure 4, the expanded range of solutions was very small when the second-order slip parameter was taken into account. In these figures, the graph of the mass transfer coefficient is not shown graphically since it completely reflected the graph of the heat transfer coefficient, because of the value of $Nb=Nt$.

Figure 5 shows the effects of the Soret ($Sr$) and Dufour (Du) parameters on the heat transfer coefficient, where the increment of the $Sr$ effect (decreased $Du$ effect) led to a decrease in the value of the heat transfer coefficient. It could also be concluded that a higher $Sr$ value decreased the heat transfer rate at the surface. Further, the increment of the $Sr$ value did not expand the range of the solutions. The effect of different values of the Brownian motion (Nb) parameter is presented in Figure 6. The heat transfer coefficient decreased with an increase in $Nb$, while the mass transfer coefficient increased with an increase in $Nb$. Therefore, the heat transfer rate at the surface increased rapidly with the smallest $Nb$, while largest values of $Nb$ were required to increase the mass transfer rate at the surface. The effects of the thermophoresis parameter ($Nt$) on the heat and mass transfer coefficients could also be seen (Figure 7). The heat transfer coefficient increased as $Nt$ increased, whereas the mass transfer coefficient decreased as $Nt$ increased. From the figures, it can be noted that the largest value of $Nt$ was required to increase the heat transfer rate, while for the mass transfer rate it was the opposite.

The velocity, temperature, and concentration profiles are presented in Figure 8, Figure 9 and Figure 10 for the effects of the slip parameters, Soret and Dufour parameters, and also for different values of $\epsilon$. It is clear from these profiles that the boundary condition (12) was converged asymptotically. Apart from that, the dual solution obtained in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 was graphically supported by the dual velocity, temperature, as well as concentration profiles presented in these figures. In addition, the boundary layer thickness for the second solution was always greater than for the first solution for each profile. The thermal and concentration boundary layer thickness increased as $Sr$ increased ($Du$ decreased). As the thermal boundary layer thickness grew larger, the heat transfer coefficient was expected to decrease (Figure 5) [59].

The stability analysis was performed to verify if either the first or second solution was stable and hence physically realizable by substituting the set of ordinary differential Equations (22)–(24) together with the respected boundary condition (25) into the codes. The purpose of the analysis was to determine the smallest eigenvalue $\gamma$ in order to identify the solutions in which a positive eigenvalue corresponded to a stable solution whereas a negative eigenvalue corresponded to an unstable solution. The first solution presented a positive value and was a stable solution characterized by a slightly disturbance in the boundary layer separation that did not interrupt the flow system (

Table 1. Smallest eigenvalues $\gamma$ for selected values of $\epsilon$ with different values of $\sigma$ and $\delta$.

σ	δ	ε	First Solution	Second Solution
0.1	−0.1	−1.257	0.0722	−0.0714
		−1.25	0.2492	−0.2397
		−1.2	0.6975	−0.6230
0.5	−0.5	−1.716	0.0295	−0.0295
		−1.71	0.1671	−0.1659
		−1.7	0.2706	−0.2676
1	−1	−2.75	0.0489	−0.0488
		−2.7	0.3561	−0.3531
		−2.6	0.6139	−0.6048

). However, the second solution presented a negative value which led to an initial growth of disturbance that interrupted the boundary layer separation and, hence, was an unstable solution. Thus, the first solution was a stable solution and, hence, it can be obtained physically because of its physical properties, whereas the second solution was an unstable solution.

6. Conclusions

The effects of the slip parameters, Soret and Dufour parameters, Brownian motion, and thermophoresis on the skin friction, heat transfer, and mass transfer coefficients in stagnation boundary-layer flow over a stretching/shrinking surface of nanofluid were investigated numerically. A dual solution was found, and the effects of the selected parameters were presented graphically. The results revealed that:

• The skin friction coefficient decreased as the first-order slip parameter and the magnitude of the second-order slip parameter ($\sigma$ and $\left|\delta \right|$) increased, whereas the heat transfer rate increased.
• The range of solutions widely expanded with the increment of the first-order slip parameter, while the expansion of the range of solutions was very small when the second-order slip parameter was considered.
• The smallest Soret number $Sr$ was required to increase the heat transfer rate at the surface.
• The heat transfer rate at the surface increased rapidly as the Brownian motion parameter $Nb$ decreased, while the largest value of $Nb$ was required to increase the mass transfer rate at the surface.
• The largest value of the thermophoresis parameter $Nt$ was required to increase the heat transfer rate, while, for the mass transfer rate, the smallest value of $Nt$ was needed.
• The first solution was a stable solution and, hence, its physical properties could be realized, whereas the second solution was an unstable solution and, hence, not physically realizable.