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

: The mathematical model of the two-dimensional steady stagnation-point ﬂow over a stretching or shrinking sheet of nanoﬂuid 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 coefﬁcient, heat transfer coefﬁcient, mass transfer coefﬁcient, as well as the velocity, temperature, and concentration proﬁles 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 (ﬁrst- 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 coefﬁcient, whereas the heat transfer coefﬁcient increased. Besides that, a larger Soret effect (smallest Dufour effect) led to the decrement of the heat transfer coefﬁcient. The effects of the Brownian motion and thermophoresis parameters on the heat transfer coefﬁcient were also studied in this paper. A stability analysis was performed in this paper to verify the stability of the solutions obtained.


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 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.

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 U w U ∞ were in the linear form, with U ∞ = 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 ∂u ∂x along with the initial and boundary conditions , v = 0, 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] where A and B are constant.

Steady-State Solution
The following similarity transformation was introduced where  is the similarity variable,  is the stream function defined as / uy     , and /      vx , which automatically satisfies Equation (1). The similarity variables (7) were substituted by Equations (2)-(4) to obtain the following ordinary (similarity) differential equations 2 10 LeSr Nb (10) subject to the boundary conditions In the above equations, primes denote the differentiation with respect to  . According to Mukhopadhyay where 0   for stretching and 0   for shrinking.
The physical quantities of practical interest are the local skin friction coefficient Nusselt number x Nu , and the local Sherwood number x Sh which are defined as:

Steady-State Solution (∂/∂t=0)
The following similarity transformation was introduced where η is the similarity variable, ψ is the stream function defined as u = ∂ψ/∂y , and v = − ∂ψ/∂x, which automatically satisfies Equation (1). The similarity variables (7) were substituted by Equations (2)-(4) to obtain the following ordinary (similarity) differential equations 1 Pr subject to the boundary conditions In the above equations, primes denote the differentiation with respect to η. According to Mukhopadhyay and Andersson [57], A = σ √ ν/c and B = δ(ν/c), with σ > 0 being the first velocity slip parameter, and δ < 0 the second velocity slip parameter (see Fang et al. [24]). Here, 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 ε is the velocity ratio parameter, which are defined as: where ε > 0 for stretching and ε < 0 for shrinking. The physical quantities of practical interest are the local skin friction coefficient C f , the local Nusselt number Nu x , and the local Sherwood number Sh x which are defined as: where τ 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 Using (7) in (13) and (14), we obtained where Re x = cx 2 /ν is the local Reynolds number.

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 τ = ct was introduced. The use of τ 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 τ and (7), we have Equations (2) and (3) can be written as subject to the boundary conditions To determine the stability of the solution, with f = f 0 (η), θ = θ 0 (η), and φ = φ 0 (η) satisfying the boundary value problem (17)-(20), we have (see Weidman et al. [44] or Roşca and Pop, [58]) where γ is an unknown eigenvalue parameter and F(η), G(η) and H(η) are small relative to f = f 0 (η), θ = θ 0 (η), and φ = φ 0 (η). The solution of the eigenvalue problem (17)- (20) gives an infinite set of eigenvalues γ 1 < γ 2 < γ 3 . . .; if γ 1 is negative, there is an initial growth of disturbances, and the flow is unstable, but when γ 1 is positive, there is an initial decay, and the flow is stable. Introducing (21) into (17)-(20), we obtain the following linearized problem subject to the boundary conditions It is inevitable to mention that for particular values of ε, σ, δ, Pr, Sr, Du, Nb, Nt, and Le, the stability of the corresponding steady flow solution f 0 (η), θ 0 (η), and φ 0 (η), was determined by the smallest eigenvalue γ. According to Harris et al. [45], the range of possible eigenvalues can be determined by relaxing a boundary condition either on F 0 (η), G 0 (η) or H 0 (η). For the present problem, we relaxed the condition that F 0 (η) → 0 , as η → ∞ and, for a fixed value of γ, we solved the system (22)-(24) subject to (25) considering the new boundary condition F 0 (0) = 1.

Results and Discussion
To discuss in detail the effects of the slip parameters (first-order slip σ and second-order slip δ), 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 (Figures 2-10). The effects of the first-order slip parameter σ and second-order slip parameter δ were plotted (Figures 2-4). From the figures, it was clear that a unique solution existed for ε > −1, dual solutions were found for ε c ≤ ε ≤ −1, and no solution was reported for ε < ε c . These three figures indicate that an increase in the slip parameters (σ and δ) 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 Figures 3 and 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  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    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 Figures 8-10 for the effects of the slip parameters, Soret and Dufour parameters, and also for different values of  . It is clear from these profiles that the boundary condition (12) was converged asymptotically. Apart from that, the dual solution obtained in Figures 2-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  in order to identify the solutions in which a positive The velocity, temperature, and concentration profiles are presented in Figures 8-10 for the effects of the slip parameters, Soret and Dufour parameters, and also for different values of ε. It is clear from these profiles that the boundary condition (12) was converged asymptotically. Apart from that, the dual solution obtained in Figures 2-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 γ 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). 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.

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 (σ and |δ|) 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.