**1. Introduction**

Nanofluids are a special class of fluids that have been the subject of developing research in the recent years. The dispersion of single nanoparticles like metals, oxides, carbon nanotubes or carbides in a fluid can create a modern class of fluids known as nanofluid. Water, oil and ethylene glycol are the common base fluid used in the formation of nanofluid. The invention of nanofluids that have good thermophysical properties can improve heat transfer performance for enormous futuristic applications such as in nuclear cooling systems, solar water heating, biomedical applications, lubrication, thermal storage, refrigeration, coolant in automobile radiator, and many others [1–13]. Choi et al. [14] initiated an experiment on nanotube-in-oil suspensions and measured that the thermal conductivity is inevitably greater than the theoretical predictions. An analytical model by Buongiorno [15] highlighted the importance of thermophoresis and Brownian motion that can induce a relative velocity between the nanoparticles and base fluid. Nield and Kuznetsov [16] implemented Buongiorno's model on the Cheng–Minkowycz problem for flow in a porous medium filled with nanofluid. A brief study of nanofluid and its thermal conductivity has also been examined by Buongiorno et al. [17]. Kuznetsov and Nield [18] revised their model [16] by introducing a new boundary condition that could manipulate the nanoparticles' volume fraction at the surface. Based on the report, it is assumed that the nanoparticles' volume fraction is passively controlled at the boundary which could make the new model more realistic and physically applicable as compared to the existing model. Muhammad et al. [19] also imposed coupled effects of convective heat and zero nanoparticles flux on conditions for Darcy–Forchheimer flow of Maxwell nanofluid. In addition, research works on the zero nanoparticles flux condition were also considered by Rehman et al. [20], Rahman et al. [21], Uddin et al. [22], ur Rahman et al. [23] and Jusoh et al. [24]. Furthermore, studies on the boundary layer problem utilizing Buongiorno's model of nanofluid were also conducted by these researchers [25–31].

Magnetohydrodynamics (MHD), also acknowledged as hydromagnetics and magneto-fluid dynamics, is the study on the behaviour of electrically conducting fluids including liquid metals, plasmas, electrolytes and salt water. The magnetic fields can generate currents or Lorentz force in a moving fluid, which give resistance to the fluid flow and, simultaneously, changes the magnetic field. Magnetohydrodynamics are effectively applied in many devices such as generators, power pumps, heat exchangers and electrostatic filters. The imposition of the magnetic field is also practical in maintaining the boundary layer flow. Rashidi et al. [32] introduced a new analytical method (DTM–Padé) to solve the boundary-layer equations of an MHD micropolar fluid near an isothermal-stretching sheet and concluded that the DTM–Padé is applicable for solving magnetohydrodynamic (MHD) boundary-layer equations. Rashidi et al. [33] concluded that the magnetic nanofluid flow over a porous rotating disk was beneficial in rotating MHD energy generators for new space systems. Sheikholeslami et al. [34] investigated the problem of an eccentric semi-annulus saturated with nanofluid under the influence of a magnetic field. Hayat et al. [35] studied the combined effects of magnetic field, velocity slip and nonlinear thermal radiation on the three-dimensional nanofluid flow. Kandasamy et al. [36] considered the convective condition for an MHD mixed convection flow in a nanofluid while Bhatti et al. [37] analyzed the MHD Wlliamson nanofluid over a porous shrinking sheet. An external magnetic field was imposed, while the induced magnetic field was neglected due to the negligible magnetic Reynolds number. Bhatti et al. [38] examined the coupled effects of MHD and partial slip on the blood flow using a Ree–Eyring fluid model. Makulati et al. [39] utilized the Tiwari and Das model of water-alumina nanofluid to study the impact of a magnetic field in inclined C-shaped enclosures. Hussain et al. [40] found that an upsurge of magnetic field, Hall current, rotation and chemical reaction had an impact on the fluid flow over an accelerated moving plate.

Stagnation point flows are fluid flows that approach the surface of a solid object and then separate into different streams. Stagnation region which meets the maximum pressure, heat transfer and mass deposition are essential in the industrial and technological field. Early classical works on stagnation flow that excluded the slip velocity effect were examined by Hiemenz [41] for the two-dimensional problem, Homann [42] for the axisymmetric case and Howarth [43] for the flow near the stagnation region. According to Wang [44], anisotropy is reflected in the difference of bi-directional (direction) slip velocities. Previously, Wang had considered the problem of two-dimensional stagnation flow with symmetric axes on dissimilar surfaces; stationary plate with isotropic slip [45] and moving plate [46], respectively. The results indicated that the enhancement of slip parameter will change the surface resistance and velocity profile. Hussain et al. [47] studied the slip flow and heat transfer of nanofluids embedded in a Darcy-type porous medium. From the analysis, the increment in the permeability of the porous medium and the velocity slip parameter increased the heat transfer rate, whereas it decreased the momentum and thermal boundary layer thicknesses. Khan et al. [48] also pointed up the significance of wall velocity slip for a reliable design and operation of microfluidic devicesmade of hydrophobic devices. The emerging numbers of industrial and technological applications make the study of anisotropic slip that relies upon the flow direction, which is influential for these types of surfaces: hydrophobic [49–54] and porous [55,56]. Rashad [57] investigated the coupled effects of anisotropic slip and convective condition in an unsteady nanofluid saturated with porous medium. Hafidzuddin et al. [58] discussed the anisotropic slip on stagnation flow due to a permeable moving surface, which resulted in dual solutions attained with the presence of a suction parameter. Numerous studies involving three-dimensional stagnation nanofluid flow towards a moving anisotropic slip surface are conducted by Raees et al. [59], Uddin et al. [60] and Balushi et al. [61]. Very recently, Sadiq [62] concluded that the intensity of magnetic field and velocity slip cause an increase in the nanofluid velocity while the boundary layer thickness decreases near the stagnation point region.

Motivated by the literature mentioned earlier, the current work accentuates the anisotropic slip effect on the numerical solution of stagnation point flow with nanoparticles due to a moving surface. Buongiorno's model of nanofluid with an assumption of zero normal flux condition at the wall is implemented. A set of transformations is applied on the governing model to simplify them into nonlinear ordinary differential equations (ODEs). The bvp4c function in MATLAB software (R2017b, MathWorks, Natick, MA, USA) is utilized to perform the numerical computations. The numerical results are demonstrated in the graph forms of velocity, nanoparticles concentration and temperature including the skin friction coefficient, heat and mass transfer rate within the specific range of related parameters. The authors also concern about the emergence of non-unique solutions and the way of stability analysis is conducted to prove the physical solution. The pioneer works on formulation of stability analysis were conducted by Merkin [63], Weidman et al. [64], Harris et al. [65] and Ro¸sca and Pop [66]. A brief of explanation on stability analysis was also discussed by the following literature [67–74]. To the best of the authors' knowledge, the results are new and have not been published.

#### **2. Mathematical Formulation**

Consider a steady, three-dimensional stagnation point flow of a nanofluid towards a moving sheet with the presence of anisotropic slip as illustrated in Figure 1. *Ux* and *Vy* are the moving plate velocities along the *x*- and *y*-directions while the *z*-axis is the axis of stagnation flow. A magnetic field of consistent strength *B*0 is applied normal to the plate. The surface is kept with fixed temperature, *Tw* while variable nanoparticles fraction, *DB ∂C∂z* + *DT T*∞ *∂T∂z* is adapted at *z* = 0. This condition is applied to achieve practically applicable results [18–24]. The ambient nanofluid concentration and temperature are denoted as *C*∞ and *T*<sup>∞</sup>.

**Figure 1.** The coordinate system of the physical model.

The following assumptions for the physical model are also contemplated in the present work:


Under all these assumptions, the flow equations are:

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0,
\tag{1}
$$

$$u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} = \mathcal{U}\_{\varepsilon}(\mathbf{x})\frac{d\mathcal{U}\_{\varepsilon}(\mathbf{x})}{dx} + \nu\nabla^2 u - \frac{\sigma\_{\mathcal{M}}\mathcal{B}\_0^{\mathcal{Q}}}{\rho} \left(u - \mathcal{U}\_{\varepsilon}(\mathbf{x})\right),\tag{2}$$

$$u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} = V\_{\varepsilon}(y)\frac{dV\_{\varepsilon}(y)}{dy} + \nu\nabla^2 v - \frac{\sigma\_M B\_0}{\rho} \left(v - V\_{\varepsilon}(y)\right),\tag{3}$$

$$
\mu \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} = \nu \nabla^2 w,\tag{4}
$$

$$u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} + w\frac{\partial T}{\partial z} = a\nabla^2 T + \tau\_1 D\_B \left[ \frac{\partial T}{\partial x}\frac{\partial C}{\partial x} + \frac{\partial T}{\partial y}\frac{\partial C}{\partial y} + \frac{\partial T}{\partial z}\frac{\partial C}{\partial z} \right] + \tau\_1 \frac{D\_T}{T\_\infty} \left[ \left(\frac{\partial T}{\partial x}\right)^2 + \left(\frac{\partial T}{\partial y}\right)^2 + \left(\frac{\partial T}{\partial z}\right)^2 \right], \tag{5}$$

$$
\mu \frac{\partial \mathcal{C}}{\partial x} + \nu \frac{\partial \mathcal{C}}{\partial y} + w \frac{\partial \mathcal{C}}{\partial z} = D\_B \nabla^2 \mathcal{C} + \frac{D\_T}{T\_{\infty}} \nabla^2 T,\tag{6}
$$

subject to the initial and boundary conditions

$$\begin{aligned} u\left(\mathbf{x},y,0\right) &= \lambda \mathcal{U}\_{\mathbf{x}} + \mu \mathcal{S}\_{1} \left. \frac{\partial u}{\partial z} \right|\_{\left(\mathbf{x},y,0\right)}, v\left(\mathbf{x},y,0\right) = \lambda V\_{\mathcal{Y}} + \mu \mathcal{S}\_{2} \left. \frac{\partial v}{\partial z} \right|\_{\left(\mathbf{x},y,0\right)}, w\left(\mathbf{x},y,0\right) = 0, \\ T\left(\mathbf{x},y,0\right) &= T\_{\mathcal{W}\_{\mathbf{y}}} \left[ D\_{\mathcal{B}} \frac{\partial \mathcal{C}}{\partial \mathbf{z}} + \frac{D\_{T}}{T\_{\infty}} \frac{\partial T}{\partial \mathbf{z}} \right]\_{\left(\mathbf{x},y,0\right)} = 0, \end{aligned} \tag{7}$$

$$\begin{aligned} \mu\left(\mathbf{x}, y, z\right) &\to \mathcal{U}\_{\varepsilon}\left(\mathbf{x}\right) = b\mathbf{x}, \quad \nu\left(\mathbf{x}, y, z\right) \to V\_{\varepsilon}\left(y\right) = by, \quad T\left(\mathbf{x}, y, z\right) \to T\_{\text{co}}\\ \mathbb{C}\left(\mathbf{x}, y, z\right) &\to \mathbb{C}\_{\infty} \quad \text{as} \quad z \to \infty, \end{aligned} \tag{8}$$

where (*<sup>u</sup>*, *v*, *w*) are the respective velocities in (*<sup>x</sup>*, *y*, *z*) directions, *T* is the fluid temperature, *C* is the nanoparticles volume fraction, *ν* is the kinematic viscosity, *μ* is the dynamic viscosity, *σM* is the electrical conductivity of the fluid, *ρ* is the fluid density, *τ*1 is the ratio of heat capacity of the nanoparticles to the base fluid, *α* is the thermal diffusivity, *DB* is the Brownian diffusion coefficient, *DT* is the thermophoretic diffusion coefficient, *S*1 is the slip coefficient in *x*-direction, *S*2 is the slip coefficient in *y*-direction, *b* is the strength of the stagnation flow, and *λ* is the moving parameter such that *λ* > 0 and *λ* < 0 refer to the moving plate, which is out and towards the origin, respectively [58,75].

The following similarity transformations which satisfy Equation (1) are employed to convert the PDEs in Equations (2)–(6) aligned with the conditions (see Equations (7) and (8)) into a set of ODEs:

$$\begin{aligned} u &= bxf'(\eta) + \mathcal{U}\_x h(\eta), \quad v = byg'(\eta) + V\_y k(\eta), \quad w = -\sqrt{b}v \left[ f(\eta) + \mathcal{g}(\eta) \right], \\ T &= \left( T\_w - T\_\infty \right) \theta(\eta) + T\_\infty, \quad \mathcal{C} = \left( \mathcal{C}\_w - \mathcal{C}\_\infty \right) \phi(\eta) + \mathcal{C}\_\infty, \quad \eta = \sqrt{\frac{b}{v}} z, \end{aligned} \tag{9}$$

where *η* is the similarity variable. Hence, the transformed nonlinear ODEs in conjunction with the conditions are:

$$f^{\prime\prime\prime} = \left(f^{\prime}\right)^2 - (f+g)f^{\prime\prime} + M\left(f^{\prime}-1\right) - 1,\tag{10}$$

$$\mathbf{g}^{\prime\prime\prime} = \left(\mathbf{g}^{\prime}\right)^2 - (f+\mathbf{g})\mathbf{g}^{\prime\prime} + M\left(\mathbf{g}^{\prime}-1\right) - 1,\tag{11}$$

$$hh'' = hf' - \left(f + \mathcal{g}\right)h' + Mh,\tag{12}$$

$$k'' = k\emptyset - \left(f + \emptyset\right)k' + Mk,\tag{13}$$

$$\boldsymbol{\theta}^{\prime\prime} = -\Pr\left[ (\boldsymbol{f} + \boldsymbol{g})\theta^{\prime} + \mathrm{N}b\theta^{\prime}\boldsymbol{\phi}^{\prime} + \mathrm{N}t(\boldsymbol{\theta}^{\prime})^2 \right],\tag{14}$$

$$\Phi^{\prime\prime} = -L\varepsilon \Pr(f+g)\phi^{\prime} - \frac{Nt}{Nb}\theta^{\prime\prime},\tag{15}$$

$$\begin{aligned} f(0) &= \mathbf{g}(0) = \mathbf{0}, \ \gamma\_1 f''(0) = f'(0), \ \gamma\_2 g''(0) = \mathbf{g}'(0), \ \lambda + \gamma\_1 h'(0) = h(0), \ \lambda + \gamma\_2 k'(0) = k(0), \\\ \theta(0) &= 1, \quad Nb\phi'(0) + Nt\theta'(0) = 0, \end{aligned} \tag{16}$$

$$f'(\eta) \to 1, \quad g'(\eta) \to 1, \quad h(\eta) \to 0, \quad k(\eta) \to 0, \quad \theta(\eta) \to 0, \quad \phi(\eta) \to 0 \quad \text{as} \quad \eta \to \infty,\\ \{17.1.1 - 1.1\} \to \{\theta(\eta) \to 0, \theta(\eta) \to 0\}, \quad \{\theta(\eta) \to 0, \theta(\eta) \to 0\}, \quad \{\theta(\eta) \to 0, \theta(\eta) \to 0\}, \quad \{\theta(\eta) \to 0, \theta(\eta) \to 0\}$$

where primes denote differentiation with respect to similarity variable *η*, *M* = *<sup>σ</sup>MB*0 *ρb* is the magnetic field parameter, *Nt* = *<sup>τ</sup>*1*DT*(*Tw* − *<sup>T</sup>*∞) *<sup>ν</sup>T*∞ is the thermophoresis parameter, *Nb* = *<sup>τ</sup>*1*DB*(*Cw* − *C* ∞) *ν* is the Brownian motion parameter, Pr = *ν α* is the Prandtl number, *Le* = *α DB* is the Lewis number, *γ*1 = *μS*1 !*b ν* and *γ*2 = *μS*2 !*b ν* are the slip parameters in the bi-directional *x*- and *y*-axis, proportionately. The physical interests in the study are the dimensionless skin friction coefficient, local Nusselt number (heat transfer rate) and local Sherwood number (mass transfer rate) which is denoted by

$$\mathbb{C}\_f Re\_x^{1/2} = f''(0), \quad \frac{Nu\_x}{Re\_x^{1/2}} = -\theta'(0), \quad \frac{Sh\_x}{Re\_x^{1/2}} = -\phi'(0), \tag{18}$$

accordingly.

#### **3. Stability Analysis**

*∂x*

The implementation of stability analysis is essential to affirm mathematically the stability and reliability of the dual solutions. The first non-unique solution which is asymptotically satisfying the boundary conditions will be denoted as the first or upper branch solution. For stability purposes, a time-dependent problem needs to be considered based on study in previous literature [67–73]. The unsteady form of Equations (2)–(6) is

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} = \mathcal{U}\_\varepsilon(\mathbf{x}) \frac{d \mathcal{U}\_\varepsilon(\mathbf{x})}{d\mathbf{x}} + \nu \nabla^2 u - \frac{\sigma\_M B\_0}{\rho} \left(u - \mathcal{U}\_\varepsilon(\mathbf{x})\right), \tag{19}$$

$$\frac{\partial v}{\partial t} + \mu \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} = V\_{\varepsilon}(y) \frac{dV\_{\varepsilon}(y)}{dy} + \nu \nabla^2 v - \frac{\sigma\_M B\_0}{\rho} \left(v - V\_{\varepsilon}(y)\right), \tag{20}$$

$$\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} = \nu \nabla^2 w,\tag{21}$$

$$\begin{split} \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} &= a \nabla^2 T + \tau\_1 D\_{\mathbb{B}} \left[ \frac{\partial T}{\partial \mathbf{x}} \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + \frac{\partial T}{\partial y} \frac{\partial \mathbf{C}}{\partial y} + \frac{\partial T}{\partial z} \frac{\partial \mathbf{C}}{\partial z} \right] \\ &+ \tau\_1 \frac{D\_T}{T\_{\infty}} \left[ \left( \frac{\partial T}{\partial \mathbf{x}} \right)^2 + \left( \frac{\partial T}{\partial y} \right)^2 + \left( \frac{\partial T}{\partial z} \right)^2 \right], \end{split} \tag{22}$$

,

$$T\_{\infty} \left[ \left< \partial \mathbf{x} \right> \right] \left< \left< \partial \mathbf{y} \right> \right> \left< \left< \partial \mathbf{z} \right> \right> \left< \left< \mathbf{\mathcal{Q}} \mathbf{z} \right> \right> \tag{23}$$
 
$$\frac{\partial \mathbf{C}}{\partial t} + u \frac{\partial \mathbf{C}}{\partial \mathbf{x}} + v \frac{\partial \mathbf{C}}{\partial \mathbf{y}} + w \frac{\partial \mathbf{C}}{\partial \mathbf{z}} = D\_{\mathbf{B}} \nabla^{2} \mathbf{C} + \frac{D\_{T}}{T\_{\infty}} \nabla^{2} T. \tag{23}$$

*∂z*

*Energies* **2019**, *12*, 1268

New transformations are applied to the unsteady problem (see Equations (19)–(23)) where *τ* is the non-dimensional time variable:

$$u = bx \frac{\partial f(\eta, \tau)}{\partial \eta} + lL\_1 h(\eta, \tau), \quad v = by \frac{\partial g(\eta, \tau)}{\partial \eta} + lyk(\eta, \tau), \quad w = -\sqrt{b\nu} \left[ f(\eta, \tau) + g(\eta, \tau) \right], \quad \left\{ \begin{aligned} &i = \sqrt{b\nu} \left[ f(\eta, \tau) + g(\eta, \tau) \right], \\ &i = \sqrt{b\nu} \left[ f(\eta, \tau) + T\_{\infty} \cdot \mathbb{C} = (\mathbb{C}\_w - \mathbb{C}\_{\infty}) \phi(\eta, \tau) + \mathbb{C}\_{\infty} \cdot \ \eta = \sqrt{\frac{b}{\nu}} z, \quad \tau = bt. \end{aligned} \right. \tag{24}$$

Using Equation (24), the following equations can be attained:

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

$$\frac{\partial^3 \mathcal{g}}{\partial \eta^3} + (f + \mathfrak{g}) \frac{\partial^2 \mathcal{g}}{\partial \eta^2} - \left(\frac{\partial \mathcal{g}}{\partial \eta}\right)^2 - M\left(\frac{\partial \mathcal{g}}{\partial \eta} - 1\right) + 1 - \frac{\partial^2 \mathcal{g}}{\partial \eta \partial \tau} = 0,\tag{26}$$

$$
\frac{
\partial^2 h
}{
\partial \eta^2
} + \left(f + \mathbf{g}\right) \frac{
\partial h
}{
\partial \eta
} - h \frac{
\partial f
}{
\partial \eta
} - Mh - \frac{
\partial h
}{
\partial \tau
} = 0,
\tag{27}
$$

$$\frac{\partial^2 k}{\partial \eta^2} + (f + \mathbf{g}) \frac{\partial k}{\partial \eta} - k \frac{\partial \mathbf{g}}{\partial \eta} - Mk - \frac{\partial k}{\partial \tau} = 0,\tag{28}$$

$$\frac{\partial^2 \theta}{\partial \eta^2} + \text{Pr} \left[ (f+g) \frac{\partial \theta}{\partial \eta} + Nb \frac{\partial \theta}{\partial \eta} \frac{\partial \phi}{\partial \eta} + Nt \left( \frac{\partial \theta}{\partial \eta} \right)^2 - \frac{\partial \theta}{\partial \tau} \right] = 0,\tag{29}$$

$$\frac{\partial^2 \phi}{\partial \eta^2} + L\varepsilon \text{Pr} \left[ (f+g)\frac{\partial \phi}{\partial \eta} - \frac{\partial \phi}{\partial \tau} \right] + \frac{Nt}{Nb} \frac{\partial^2 \theta}{\partial \eta^2} = 0,\tag{30}$$

with the transformed conditions

$$\begin{aligned} f(0,\tau) &= g(0,\tau) = 0, \quad \left[\gamma\_1 \frac{\partial^2 f}{\partial \eta^2} - \frac{\partial f}{\partial \eta}\right]\_{(0,\tau)} = 0, \quad \left[\gamma\_2 \frac{\partial^2 g}{\partial \eta^2} - \frac{\partial g}{\partial \eta}\right]\_{(0,\tau)} = 0, \\\ \lambda + \gamma\_1 \frac{\partial h}{\partial \eta}\bigg|\_{(0,\tau)} &= h(0,\tau), \quad \lambda + \gamma\_2 \frac{\partial k}{\partial \eta}\bigg|\_{(0,\tau)} = k(0,\tau), \quad \theta(0,\tau) = 1, \quad \left[Nb\frac{\partial \phi}{\partial \eta} + Nt\frac{\partial \theta}{\partial \eta}\bigg|\_{(0,\tau)} = 0,\right] \end{aligned} \tag{31}$$

$$\left. \frac{\partial f}{\partial \eta} \right|\_{\left(\eta,\tau\right)} \to 1, \left. \frac{\partial g}{\partial \eta} \right|\_{\left(\eta,\tau\right)} \to 1, \left. h(\eta,\tau) \to 0, \left. k(\eta,\tau) \to 0, \left. \theta(\eta,\tau) \to 0, \left. \phi(\eta,\tau) \to 0 \text{ as } \eta \to \infty \right. \right. \right| \tag{32}$$

For the stability process, steady flow solutions *f*(*η*) = *f*0(*η*), *g*(*η*) = *g*0(*η*), *h*(*η*) = *h*0(*η*), *k*(*η*) = *k*0(*η*), *<sup>θ</sup>*(*η*) = *<sup>θ</sup>*0(*η*) and *φ*(*η*) = *φ*0(*η*) which have satisfied Equations (10)–(17) are examined by the following expressions:

$$\begin{cases} f(\eta,\tau) = f\_0(\eta) + e^{-\sigma\tau} F(\eta) \\ g(\eta,\tau) = g\_0(\eta) + e^{-\sigma\tau} G(\eta) \\ h(\eta,\tau) = h\_0(\eta) + e^{-\sigma\tau} H(\eta) \\ k(\eta,\tau) = k\_0(\eta) + e^{-\sigma\tau} K(\eta) \\ \theta(\eta,\tau) = \theta\_0(\eta) + e^{-\sigma\tau} P(\eta) \\ \phi(\eta,\tau) = \phi\_0(\eta) + e^{-\sigma\tau} R(\eta) \end{cases} \tag{33}$$

where *σ* is an eigenvalue, *<sup>F</sup>*(*η*), *<sup>G</sup>*(*η*), *<sup>H</sup>*(*η*), *<sup>K</sup>*(*η*), *<sup>P</sup>*(*η*) and *<sup>R</sup>*(*η*) are small relative to *f*0(*η*), *g*0(*η*), *h*0(*η*), *k*0(*η*), *<sup>θ</sup>*0(*η*) and *φ*0(*η*), correspondingly [66]. The linearized eigenvalue problem will be attained by replacing Equation (33) into Equations (25)–(32):

$$F'''' + \left(f\_0 + g\_0\right)F'' + \left(F + G\right)f\_0'' - \left(2f\_0' - \sigma + M\right)F' = 0,\tag{34}$$

$$\left(\mathbf{G}''' + \left(f\_0 + \mathbf{g}\_0\right)\mathbf{G}'' + \left(\mathbf{F} + \mathbf{G}\right)\mathbf{g}\_0'\right) - \left(2\mathbf{g}\_0' - \sigma + M\right)\mathbf{G}' = 0,\tag{35}$$

$$H'' + \left(f\_0 + \mathfrak{g}\_0\right)H' + \left(F + G\right)h\_0' - h\_0F' - \left(f\_0' - \sigma + M\right)H = 0,\tag{36}$$

$$\mathcal{K}^{\prime\prime} + \left(f\_0 + \mathcal{g}\_0\right)\mathcal{K}^{\prime} + \left(F + G\right)k\_0^{\prime} - k\_0\mathcal{G}^{\prime} - \left(\mathcal{g}\_0^{\prime} - \sigma + M\right)\mathcal{K} = 0,\tag{37}$$

$$\Pr^{\prime\prime} + \Pr\left[\left(f\_0 + \mathfrak{g}\_0\right)P^{\prime} + \left(F + G\right)\theta\_0^{\prime} + \mathrm{Nb}\left(\mathfrak{gl}\_0^{\prime}P^{\prime} + \theta\_0^{\prime}R^{\prime}\right) + \mathrm{Nt}\left(2\theta\_0^{\prime}P^{\prime}\right) + \sigma P\right] = 0,\tag{38}$$

$$\mathcal{R}^{\prime\prime} + L\boldsymbol{x}\Pr\left[\left(f\mathrm{o} + \mathrm{go}\right)\mathcal{R}^{\prime} + \left(\mathrm{F} + \mathrm{G}\right)\phi\_{0}^{\prime} + \sigma\mathcal{R}\right] + \frac{N\hbar}{N\hbar}\mathcal{P}^{\prime\prime} = 0,\tag{39}$$

along with the conditions

$$\begin{aligned} F(0) &= G(0) = 0, & \gamma\_1 F''(0) &= F'(0), & \gamma\_2 G''(0) &= G'(0), \\ \gamma\_1 H'(0) &= H(0), & \gamma\_2 K'(0) &= K(0), & P(0) &= 0, & NbR'(0) + NtP'(0) &= 0 \end{aligned} \tag{40}$$

$$F'(\eta) \to 0, \ G'(\eta) \to 0, \ H(\eta) \to 0, \ K(\eta) \to 0, \ P(\eta) \to 0, \ R(\eta) \to 0 \qquad \text{as} \qquad \eta \to \infty. \tag{41}$$

The stability of the solutions *f*0(*η*), *g*0(*η*), *h*0(*η*), *k*0(*η*), *<sup>θ</sup>*0(*η*) and *φ*0(*η*) depends on the smallest eigenvalue, *σ*1 by solving the governing linearized eigenvalue model in Equations (34)–(41). Relaxation of a boundary condition is necessary to attain a possible range of eigenvalues [65]. Hence, in the present paper, the boundary condition *<sup>F</sup>*(*η*) → 0 as *η* → ∞ (see Equation (41)) is relaxed and replaced with the normalizing boundary condition *F*(0) = 1.

#### **4. Results and Discussion**

The similarity solutions for the transformed ODEs in Equations (10)–(15) aligned with the conditions (see Equations (16) and (17)) are found with the aid of bvp4c function in MATLAB software. The bvp4c function implements a finite difference scheme known as 3-stage Lobatto IIIa [74,76–78]. There are four separate codes in bvp4c function; code a for solving steady flow equations, code b for continuation of code a, code c and d for stability analysis of dual solutions. In the bvp4c code a, *η*∞ = 10 is used for the numerical calculations, but it is found that *η*∞ = 5 is sufficient enough to fulfill the boundary conditions based on all the profiles demonstrated in the present study. For the method validation, few values in the case of Newtonian fluid with the absence of magnetic parameter and heat transfer have been compared to the results by Wang [44] and Balushi et al. [61]. The present numerical values are in positive agreemen<sup>t</sup> with others as tabulated in Table 1.


**Table 1.** A comparison data with previous published results for *λ* = 1, *γ*1 = 5 and various *γ*2.

Figures 2–17 exhibit that dual solutions are possible in the present study. The first and second solutions are depicted by the straight and dashed lines, respectively. Generally, the solution which converges first is assumed as the first or physical solution, and in this case, the stability analysis as discussed in the previous section will affirm which solution is physically realizable. Figures 2–7 display the nanofluid velocities on *x*- and *y*-directions, temperature and concentration profiles in limiting numbers of *γ*2 when *Le* = Pr = 1, *Nb* = *Nt* = 0.2, *M* = 0.5, *λ* = −1, 1 and *γ*1 = 5. Both of the velocities as revealed in Figures 2 and 3 increase because the shear stress decrease with the enlargement of the slip parameter and this is supported by the result in Figure 11. Figures 4 and 5 demonstrate the slip velocity profiles *h*(*η*) and *k*(*η*) in both directions, respectively. There are opposite and symmetrically effect for both profiles at *λ* = ±1. The changes in *λ* only influence the slip velocity profiles since the moving parameter only remains in the initial condition at *h*(0) and *k*(0) (see Equation (16)). Boundary layer thickness for both solutions diminishes with an increment of *γ*2 as manifested in Figures 4 and 5.

The nanofluid temperature as illustrated in Figure 6 declines with the expanding values of *γ*2 while the opposite result is obtained for the nanoparticles' volume fraction profile in Figure 7. Meanwhile, Figures 8–10 present the influence of the *Nb* and *Nt* on both nanoparticles volume fraction (concentration) and temperature profiles. The nanofluid concentration diminishes as *Nb* increases while expands as *Nt* increases and these outcomes are in conjunction with the previous results by Balushi et al. [61]. An inflation of *Nt* seems to boost the thermal boundary layer thickness and temperature, whereas the presence of *Nb* give zero impact on the nanofluid temperature. The appearance of the nanoparticles in the base fluid will generate thermophoresis parameter. As *Nt* increases, the higher thermal conductivity in nanofluid will enhance both temperature and concentration.

Variations of *g*(0), *<sup>k</sup>*(0), −*θ*(0) and −*φ*(0) against the slip parameter, *γ*1 for selected values of *γ*2, *Nb* and *Nt* are visualized in Figures 11–17, respectively. Unlike the other studies which conduct stability analysis for non-unique solutions, no critical value or turning point is found in this analysis. Critical value is defined as a value that separates the first (physical) and second solutions. Both upper and lower branches existing boundlessly for each value of *γ*1 and for the graph visualization, 0 ≤ *γ*1 ≤ 30 have been selected. There is no significant effect of *γ*2 on *f* (0) and *<sup>h</sup>*(0), hence the graph is not highlighted. This is supported by the velocity and slip velocity profile that have been shown previously in Figures 2 and 4. Variations of −*θ*(0) and −*φ*(0) as can be seen in Figures 13 and 14 only show minimal changes as *γ*2 increase but since the value of *Nb* and *Nt* is same, −*θ*(0) = −(−*φ*(0)) for all values of *γ*1 and *γ*2. An increasing values of the Brownian motion parameter *Nb* escalate the mass transfer rate whereas the thermophoresis parameter *Nt* reduces both heat and mass transfer rate as displayed in Figures 15–17.

Since dual solutions are obtained in the study, a stability analysis has been conducted by solving Equations (34)–(39) with the conditions (see Equations (40) and (41)) using Matlab bvp4c code c and d. Generally, negative *σ*1 indicates an initial upsurge of disturbances, which signifies that the flow is unstable while there are opposite flow characteristics for positive *σ*1. The smallest eigenvalue, *σ*1 for some values of *γ*1, is tabulated in Table 2. It shows that the first and second solutions have positive and negative eigenvalues, respectively, which validate the stability and reliability of the first solution.

**Figure 2.** Nanofluid velocity along the *x*-direction for diverse values of *γ*2 with *λ* = −1, 1.

**Figure 3.** Nanofluid velocity along the *y*-direction for diverse values of *γ*2 with *λ* = −1, 1.

**Figure 4.** Slip velocity along the *x*-direction for diverse values of *γ*2 with *λ* = −1, 1.

**Figure 5.** Slip velocity along the *y*-direction for diverse values of *γ*2 with *λ* = −1, 1.

**Figure 6.** Nanofluid temperature in a limiting range of *γ*2 with *λ* = −1, 1.

**Figure 7.** Nanoparticles volume fraction profile in a limiting range of *γ*2 with *λ* = −1, 1.

**Figure 8.** Nanoparticles' volume fraction profile in limiting range of *Nb* when *Le* = Pr = 1, *Nt* = 0.2, *γ*1 = 5, *γ*2 = 10, *M* = 0.5 and *λ* = −1, 1.

**Figure 9.** Nanofluid temperature in limiting range of *Nt* when *Le* = Pr = 1, *Nb* = 0.2, *γ*1 = 5, *γ*2 = 10, *M* = 0.5 and *λ* = −1, 1.

**Figure 10.** Nanoparticles' volume fraction profile in a limiting range of *Nt* when *Le* = Pr = 1, *Nb* = 0.2, *γ*1 = 5, *γ*2 = 10, *M* = 0.5 and *λ* = −1, 1.

**Figure 11.** Variations of *g*(0) towards *γ*1 for selected values of *γ*2.

**Figure 12.** Variations of *k*(0) towards *γ*1 for selected values of *γ*2.

**Figure 13.** Variations of −*θ*(0) towards *γ*1 for selected values of *γ*2.

**Figure 14.** Variations of −*φ*(0) towards *γ*1 for selected values of *γ*2.

**Figure 15.** Variations of −*φ*(0) towards *γ*1 for selected values of *Nb*.

**Figure 16.** Variations of −*θ*(0) towards *γ*1 for selected values of *Nt*.

 **Figure 17.** Variations of −*φ*(0) towards *γ*1 for selected values of *Nt*.


**Table 2.** *σ*1 for diverse *γ*1 when *γ*2 = 5, *Le* = 1, *Nb* = *Nt* = 0.2, *M* = 0.5, Pr = 1 and *λ* = −1.

## **5. Conclusions**

A numerical investigation on three-dimensional MHD stagnation point flow due to a moving plate with the presence of anisotropic slip has been accomplished. Buongiorno's model of nanofluid is selected to integrate the thermophoresis and Brownian motion parameters. A set of transformations was used to reduce the governing model into a set of nonlinear differential equations. The similarity equations are then transformed into the bvp4c algorithm to perform the numerical computation. The numerical results are illustrated graphically on the specific physical parameters such as thermophoresis parameter *Nt*, slip parameters *γ*1, *γ*2 and Brownian motion parameter *Nb* to study the flow, heat and mass transfer characteristics of the nanofluid. Non-unique solutions exist boundlessly for all positive values of related parameters in the present study. The implementation of stability analysis using bvp4c code proves the reliability of the first solution. The thickness of boundary layer for the non-physical solution is larger compared to the physical solution. The rate of heat transfer is reduced with the augmentation of the anisotropic slip (difference of slip velocities), while the opposite result is obtained for the mass transfer rate. The rising values of *Nb* reduce the boundary layer thickness for nanoparticle concentration and increase the Sherwood number while an adverse trend is observed for the higher values of *Nt*.

**Author Contributions:** Research design, N.S.K., N.M.A., R.N. and I.P.; Formulation and methodology, N.S.K. and N.M.A.; Result analysis, N.S.K.; Validation, N.M.A., R.N., E.H.H. and N.W.; Article preparation, N.S.K.; Review and editing, N.M.A., R.N., E.H.H., N.W. and I.P.

**Funding:** The authors would like to express their grea<sup>t</sup> appreciation to the Universiti Putra Malaysia through the Putra Grant (9570600). The main author also would like to acknowledge Universiti Teknikal Malaysia Melaka and Ministry of Education (Malaysia) for the financial support through a UTEM-SLAB scholarship. The work by Ioan Pop has been supported from the gran<sup>t</sup> PN-III-P4-ID-PCE-2016-0036, UEFISCDI, Romania.

**Acknowledgments:** The authors appreciate the valuable feedback and recommendations by the reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest. *Energies* **2019**, *12*, 1268
