**1. Introduction**

Mixed convection flows or a combination of forced and free convections exists in numerous transport operations, both naturally occurring and in engineering applications. Such applications for example, include heat exchangers, solar collectors, nuclear reactors, atmospheric boundary layer flow, nanotechnology, electronic apparatus, etc. These operations occur during the effects of buoyancy forces in forced convections or the effects of forced flow in free convections become substantial. Over the past several decades, most research in mixed convection flow analysis has emphasised the occurrence of dual solutions for a particular range of the buoyancy (mixed convection) parameter in the opposing flow region, such as in the research by Ramachandran et al. [1], Merkin and Mahmood [2], Devi et al. [3] and Lok et al. [4]. In contrast to [1–4], Ridha and Curie [5] continued the study by Merkin and Mahmood [2] by establishing the existence of dual solutions in both the opposing and assisting flow regions. Furthermore, by implementing a stability analysis of the dual solutions for

mixed convection flow in a saturated porous medium, Merkin [6] demonstrated that the upper branch of the solutions is stable whereas the lower branch shows instability. Accordingly, various other researches have also stated the occurrence of dual solutions in the mixed convection flow in different configurations, namely by Ro¸sca et al. [7], Rahman et al. [8] and recently by Abbasbandy et al. [9]. An inclusive account of the theoretical research prior to 1987 for both laminar and turbulent mixed convection boundary layer flows may be found in the books by Gebhart et al. [10], Schlichting and Gersten [11], Pop and Ingham [12] and Bejan [13], for example.

The innovative idea of nanofluids was first brought up by Choi et al. [14] in 1995, when the authors suggested a path for exceeding the performance of heat transfer fluids which were currently available. An extraordinary enhancement in the thermal properties of base fluids may be achieved just by utilizing a minimal amount of nanoparticles scattered uniformly and suspended stably in a base fluid. Nanofluids, as colloidal mixtures of nanoparticles (1–100 nm) along with a base liquid (nanoparticle fluid suspensions) are known, provide access to a new class of nanotechnology-based heat transfer media (Das et al. [15]). Numerous techniques and methodologies, such as rising either the heat transfer surface or the heat transfer coefficient between the fluid and the surface that allows high heat transfer rates in a small volume, may be utilized to promote heat transfer. Notwithstanding, cooling turns out to be one of the most critical technical challenges faced by numerous and diverse industries, including microelectronics, transportation, solid-state lighting, and manufacturing. The addition of micrometreor millimetre-sized solid metal or metal oxide particles to base fluids produces an increase in the thermal conductivity of the resultant fluids. On the other hand, apart from being applied in the field of heat transfer, nanofluids (nanometre-sized particles in a fluid) may also be synthesised for unique magnetic, electrical, chemical, and biological applications (see Manca et al. [16]). Nanoparticles are produced from various materials such as copper (Cu), alumina (Al2O3), titania (TiO2), copper oxide (CuO) as well as silver (Ag) (see Oztop and Abu-Nada [17]). References on nanofluids are mentioned in the books written by Das et al. [15], Nield and Bejan [18], Minkowycz et al. [19] and Shenoy et al. [20], and also in the review papers written by Buongiorno et al. [21], Kakaç and Pramuanjaroenkij [22], Fan and Wang [23], Mahian et al. [24], Sheikholeslami and Ganji [25], Gro¸san et al. [26], Myers et al. [27], etc. These review papers elaborate specifically on the production of nanofluids, the theoretical and experimental exploration of the thermal conductivity and viscosity of nanofluids, as well as the work conducted on the convective transport of nanofluids.

Interestingly, many studies investigating the boundary layer problem of mixed convection flow in a nanofluid are reported in the literature. Tamim et al. [28] examined the effects of the magnetic field, suction/injection and solid volume fraction of nanoparticles on mixed convection about the stagnation-point flow of a nanofluid. On the other hand, Subhashini et al. [29] investigated the mixed convection flow about the stagnation-point region over an exponentially stretching/shrinking sheet in a nanofluid for both suction and injection cases. Later, Mustafa et al. [30] extended the study conducted by Tamim et al. [28] in consideration of the combined effects of viscous dissipation and the magnetic field by gaining a unique solution for assisting and opposing flow cases. Recently, Ibrahim et al. [31], Mabood et al. [32] and Othman et al. [33], similarly investigated the problem of mixed convection boundary layer flow in nanofluids under different physical conditions.

The impact of thermal radiation on heat transfer becomes increasingly important in the design of advanced energy conversion systems operating at high temperature. Moreover, thermal radiation has applications in numerous technological problems such as combustion, nuclear reactor safety, solar collectors, furnace design and many others (see Ozisik [34]). Furthermore, the study of thermal radiation on flow and heat transfer characteristics in a nanofluid have attracted immense interest because nanofluids have different properties than those found in either the particles or the base fluid. Given this fact, many researchers have explored the impact of thermal radiation on flow and heat transfer in a nanofluid along with other various aspects. An important analysis by Hady et al. [35] studied the boundary layer viscous flow and heat transfer characteristics of a nanofluid over a nonlinearly stretching sheet in the presence of thermal radiation in a single-phase model. In a separate

study, Ibrahim and Shankar [36] investigated the influences of thermal radiation, magnetic fields and slip boundary conditions on boundary layer flow and heat transfer past a permeable stretching sheet in a nanofluid. Notwithstanding, Haq et al. [37] discussed the combined effects of thermal radiation, magnetohydrodynamic (MHD), velocity and thermal slip on the boundary layer stagnation-point flow of nanofluid and the effects over a stretching sheet. More recently, Daniel et al. [38] investigated the effects of thermal radiation, magnetic fields, electrical fields, Ohmic dissipation, thermal and concentration stratifications on the flow and heat transfer of electrically conducting nanofluid past a permeable stretching sheet. In another recent study, Sreedevi et al. [39] analysed the effect of thermal radiation, magnetic field and the chemical reaction on flow, heat and mass transfer of nanofluid over a linear and nonlinear stretching sheet saturated by the porous medium. Accordingly, several other studies have been undertaken on mixed convection boundary layer flow in nanofluids in the presence of thermal radiation, including works by Yazdi et al. [40], Pal and Mandal [41] and Ayub et al. [42].

The heat source/sink effect in addition to the thermal radiation effect plays a vital role in governing the heat transfer in industrial operations in which the attributes of the output are dependent on the factors of heat control. Accordingly, many researchers have studied the impacts of a heat source/sink on the boundary layer flow and heat transfer of nanofluids along with different aspects. Rana and Bhargava [43] numerically investigated the impact of the various types of nanoparticles on mixed convection flow of nanofluid along the vertical plate with a heat source/sink. Furthermore, Pal et al. [44] analysed the combined impacts of internal heat generation/absorption, thermal radiation and suction/injection on mixed convection stagnation point flow of nanofluids over a stretching/shrinking sheet in a porous medium. In addition, Pal and Mandal [45] discussed the impacts of microrotation and nanoparticles on boundary layer flow in nanofluids in the occurrence of non-uniform heat source/sink, suction, thermal radiation and magnetic fields. In another paper, Mondal et al. [46] considered the influence of heat generation/absorption and thermal radiation on hydromagnetic three-dimensional mixed convection flow of nanofluid over a vertical stretching surface. Sharma and Gupta [47] further investigated the effect of heat generation/absorption, MHD, thermal radiation, viscous dissipation on flow and heat transfer of Jeffrey nanofluids.

Interestingly, previous studies did not include the combined effects of thermal radiation, heat source/sink and suction on mixed convection flow of a nanofluid. Therefore, the primary aim of this article is to examine the impact of thermal radiation, heat source/sink and suction on mixed convection stagnation point flow over a stretching/shrinking sheet in a nanofluid, by applying a mathematical nanofluid model suggested by Tiwari and Das [48]. In our opinion, the problem is relatively new, novel with no such articles reported at this stage in the literature. Suitable similarity transformations are employed to transform nonlinear partial differential equations into nonlinear ordinary differential equations. The equations are then solved numerically with the assistance of the bvp4c programme in MATLAB, and the results are graphically plotted and displayed in tables. The results from the study indicates that dual solutions exist for a particular range of parameters, namely, the mixed convection parameter, solid volume fraction of nanoparticles, thermal radiation parameter, heat source/sink parameter, suction parameter and stretching/shrinking parameter. Further, it is useful to mention, that the stability analysis of the dual solutions is conducted to investigate which solution is stable.

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

This study considers the two-dimensional steady mixed convection flow of a viscous and incompressible nanofluid near the stagnation-point past a permeable vertical stretching/shrinking surface with the velocity *uw*(*x*) and free stream velocity *ue*(*x*), as illustrated in Figure 1, where *x* and *y* denote the Cartesian coordinates evaluated along the surface of the stretching/shrinking sheet and normal to it, respectively. The fluid consists of a water-based nanofluid comprising three distinct types of nanoparticles which are copper (Cu), alumina (Al2O3) and titania (TiO2). The thermophysical properties of water (the base fluid) and nanoparticles are shown in Table 1. These thermophysical properties will be used in the numerical computations of this study. Also, it is assumed that the

flow was subjected to the combined impact of thermal radiation and a heat source/sink. Another assumption made are such that the temperature of the stretching/shrinking sheet, *Tw*(*x*), and the temperature of the ambient nanofluid adopt a constant value *T*<sup>∞</sup>. Furthermore, it is also assumed that the water-based fluid and the nanoparticles are in thermal equilibrium and that no slip exists among them. The mathematical nanofluid model suggested by Tiwari and Das [48] is applied in this case. It should be mentioned that this nanofluid model is a single-phase approach where the nanoparticles are assumed to have a uniform shape and size, and the interactions between nanoparticles and surrounding fluid are also neglected (Pang et al. [49], Ebrahimi et al. [50] and Sheremet et al. [51]). This assumption is practical when the base fluid is easily fluidized, so it can be considered to behave as a single fluid, hence it applies to the justification of using single phase model in this study.

**Figure 1.** Physical model and coordinate system. (**a**) Stretching surface; (**b**) Shrinking surface.


**Table 1.** Thermophysical properties of water and nanoparticles (Oztop and Abu-Nada [17]).

By taking into considerations of these assumptions together with the Boussinesq and the boundary layer approximations, the governing boundary layer equations of continuity, momentum and thermal energy in the existence of thermal radiation and the heat source or sink, are given as shown below:

$$
\frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} = 0,
\tag{1}
$$

$$u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = u\_{\varepsilon}\frac{du\_{\varepsilon}}{dx} + \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2 u}{\partial y^2} + \frac{(\rho\beta)\_{nf}}{\rho\_{nf}}(T - T\_{\infty})\mathbf{g}\_{\prime} \tag{2}$$

$$u\frac{\partial T}{\partial \mathbf{x}} + v\frac{\partial T}{\partial y} = a\_{nf}\frac{\partial^2 T}{\partial y^2} - \frac{1}{\left(\rho \mathbb{C}\_p\right)\_{nf}}\frac{\partial q\_r}{\partial y} + \frac{Q\_0}{\left(\rho \mathbb{C}\_p\right)\_{nf}}(T - T\_\infty),\tag{3}$$

and the associated boundary conditions to present the flow are:

$$\begin{aligned} u &= u\_{\mathfrak{w}}(\mathbf{x}), \; v = v\_{\mathfrak{w}}(\mathbf{x}), \; T = T\_{\mathfrak{w}}(\mathbf{x}) \; \text{at } y = 0, \\ u &\to u\_{\mathfrak{e}}(\mathbf{x}), \; T \to T\_{\mathfrak{e}\mathfrak{e}}(\mathbf{x}) \; \text{as } y \to \mathfrak{so}, \end{aligned} \tag{4}$$

where *u* and *v* represent the velocity elements along the *x* and *y* directions, respectively; *g* stands for the acceleration caused by gravity, *T* denotes the temperature of the nanofluid, *Q*0 denotes the heat source/sink coefficient, with *Q*0 > 0 corresponding to the heat source and *Q*0 < 0 corresponding to the heat sink. Further, *vw*(*x*) represents the wall mass flux, with *vw*(*x*) < 0 corresponding to the suction. Moreover, *μn f* represents the dynamic viscosity of the nanofluid, *ρn f* refers to the density of the nanofluid, (*ρβ*)*n f* denotes the thermal expansion coefficient of the nanofluid as described in the Brinkman's model, *<sup>ρ</sup>Cpn f* denotes the heat capacitance of the nanofluid, *αn f* denotes the thermal diffusivity of the nanofluid, *qr* represents the radiation heat flux, and lastly, *νn f* reflects the kinematic viscosity of the nanofluid. The relations of *μn f* , *αn f* , *ρn f* , (*ρβ*)*n f* , *<sup>ρ</sup>Cpn f* and *kn f* are described in the following equations (see Oztop and Abu-Nada [17]):

$$\begin{aligned} \mu\_{nf} &= \frac{\mu\_{nf}}{(1-\phi)^{2.5}}, \; a\_{nf} = \frac{k\_{nf}}{\left(\rho \mathbb{C}\_{p}\right)\_{nf}}, \; \rho\_{nf} = (1-\phi)\rho\_{f} + \phi\rho\_{s}, \\ \left(\rho\beta\right)\_{nf} &= (1-\phi)\left(\rho\beta\right)\_{f} + \phi\left(\rho\beta\right)\_{s'}\left(\rho\mathbb{C}\_{p}\right)\_{nf} = (1-\phi)\left(\rho\mathbb{C}\_{p}\right)\_{f} + \phi\left(\rho\mathbb{C}\_{p}\right)\_{s'} \end{aligned} \tag{5}$$
 
$$\frac{k\_{nf}}{k\_{f}} = \frac{k\_{s} + 2k\_{f} - 2\phi\left(k\_{f} - k\_{s}\right)}{k\_{s} + 2k\_{f} + \phi\left(k\_{f} - k\_{s}\right)},$$

where *μf* refers to the dynamic viscosity of the base fluid, *φ* denotes the solid volume fraction of the nanoparticles, *ρf* and *ρs* represent the density of the base fluid and the density of the solid nanoparticle, respectively, *kn f* represents the thermal conductivity of the nanofluid, as approximated by the Maxwell-Garnett's model, the subscript '*f*' represents the base fluid, and lastly, '*s*' reflects the solid nanoparticle.

Meanwhile, upon employing the Rosseland's approximation, the radiation heat flux, *qr* is given by Zheng [52], which adopts the following form:

$$q\_{l'} = -\frac{4\sigma^\*}{3k^\*} \frac{\partial T^4}{\partial y} \,, \tag{6}$$

where *σ*<sup>∗</sup> represents the Stefan-Boltzmann constant and *k*∗ is the Rosseland mean spectral absorption coefficient. Furthermore, it is assumed that the temperature difference between the flow is such that *T*<sup>4</sup> can be expanded using Taylor's series as a linear combination of the temperature. Next, after the expansion of *T*<sup>4</sup> into the Taylor's series for *T*<sup>∞</sup>, the approximation was obtained by omitting the higher order terms, obtaining *T*<sup>4</sup> = <sup>4</sup>*T*3∞*<sup>T</sup>* − <sup>3</sup>*T*4∞. Therefore, upon substituting Equations (5) and (6) into Equation (3), the following equation is obtained:

$$
\mu \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial y} = a\_{nf} \frac{\partial^2 T}{\partial y^2} + \frac{16 \nu^\* T\_{\infty}^3}{3k^\* \left(\rho \mathbf{C}\_p\right)\_{nf}} \frac{\partial^2 T}{\partial y^2} + \frac{Q\_0}{\left(\rho \mathbf{C}\_p\right)\_{nf}} (T - T\_{\infty}). \tag{7}
$$

To determine the similar forms of the Equations (1), (2) and (7), with boundary conditions (4), the terms are defined; *uw*(*x*), *vw*(*x*), *Tw*(*x*) and *ue*(*x*) in the following form:

$$\mathbf{u}\_{\text{w}}(\mathbf{x}) = b\mathbf{x}, \ v\_{\text{w}}(\mathbf{x}) = -\sqrt{\mathbf{a}\mathbf{v}\_{f}}\mathbf{s}, \ T\_{\text{w}}(\mathbf{x}) = T\_{\text{w}} + T\_{0}\mathbf{x}, \ u\_{\text{\textdegree}}(\mathbf{x}) = a\mathbf{x}. \tag{8}$$

Here, *a* and *b* are constants, *s* denotes the suction parameter and *T*0 represents the constant characteristic temperature, with *T*0 < 0 indicating the cooled surface (opposing flow) while *T*0 > 0 signifies the heated surface (assisting flow).

Furthermore, the governing Equations (1), (2) and (7) together with the boundary conditions (4) have been transformed into ordinary differential equations by the dimensionless functions *u*, *v* and *θ*, in relation to the suitable similarity variable *η* as follows:

$$u = a x f'(\eta), \ v = -\sqrt{a \nu\_f} f(\eta), \ \theta(\eta) = \frac{T - T\_{\infty}}{T\_w - T\_{\infty}}, \ \eta = \sqrt{\frac{a}{\nu\_f}} y. \tag{9}$$

Note that *f*(*η*) denotes the dimensionless stream function, *f* (*η*) be the dimensionless velocity profile, *<sup>θ</sup>*(*η*) represents the dimensionless temperature profile and the prime indicates the differentiation with respect to *η*.

Equation (1) is therefore satisfied identically with the given similarity transformation (9). After substituting similarity transformation (9) into Equations (2) and (7), we obtain the following coupled nonlinear ordinary differential equations:

$$\frac{1}{\left(1-\phi\right)^{2.5}}f''' + \left(1-\phi+\phi\frac{\rho\_s}{\rho\_f}\right)\left(f f'' + 1 - \left(f'\right)^2\right) + \left(1-\phi+\phi\frac{\left(\rho\beta\right)\_\varepsilon}{\left(\rho\beta\right)\_f}\right)\lambda\theta = 0,\tag{10}$$

$$\frac{1}{\Pr} \left( \frac{k\_{nf}}{k\_f} + \frac{4}{3} \text{Nr} \right) \theta'' + \left( 1 - \phi + \phi \frac{\left( \rho \mathbb{C}\_p \right)\_s}{\left( \rho \mathbb{C}\_p \right)\_f} \right) \left( f \theta' - f' \theta + K \theta \right) = 0,\tag{11}$$

while the boundary conditions (4) adopt the new form:

$$\begin{array}{c} f(\eta) = \text{s, } f'(\eta) = \text{c, } \theta(\eta) = 1 \text{ at } \eta = 0, \\\ f'(\eta) = 1, \; \theta(\eta) = 0 \text{ as } \eta \to \infty, \end{array} \tag{12}$$

where Pr denotes the Prandtl number, *λ* represents the mixed convection parameter with the case of *λ* < 0 corresponds to the opposing flow, whereas *λ* > 0 corresponds to the assisting flow. Moreover, *Nr* denotes the thermal radiation parameter, *K* represents the heat source/sink parameter with the case *K* > 0 refers to the heat source and *K* < 0 refers to the heat sink. Further *c* denotes the stretching/shrinking parameter, with *c* > 0 for a stretching sheet and *c* < 0 for a shrinking sheet, and *s* is the constant mass flux parameter, with *s* > 0 for suction and *s* < 0 for injection or withdrawal of the fluid, The parameters *Pr*, *λ*, *Nr*, *K*, *s* and *c* can be expressed in the following equations as:

$$\Pr = \frac{\nu\_f}{a\_f}, \lambda = \frac{Gr\_\mathbf{x}}{\mathrm{Re}\_\mathbf{x}^2}, \operatorname{Nr} = \frac{4T\_\infty^3 \sigma^\*}{k\_f k^\*}, \operatorname{K} = \frac{\operatorname{Q}\_0}{a \left(\rho \mathbb{C}\_p\right)\_{\mathrm{nf}}}, \operatorname{s} = -\frac{\nu\_\mathrm{w} \left(\mathbf{x}\right)}{\sqrt{\pi \nu\_f}}, \operatorname{c} = \frac{\mathrm{b}}{a}. \tag{13}$$

Here, the local Grashof number *Grx* and the local Reynolds number Re*x* are given by:

$$\mathrm{Gr}\_x = \frac{\mathrm{g}\beta\_f (T\_w - T\_\infty) x^3}{\nu\_f^2}, \ \mathrm{Re}\_x = \frac{\mu\_\mathrm{c} x}{\nu\_f}. \tag{14}$$

The interested physical quantities are the skin friction coefficient *Cf* and the local Nusselt number *Nux* which are expressed by:

$$\mathcal{C}\_{f} = \frac{\mathfrak{T}\_{w}}{\rho\_{f}\mathfrak{u}\_{\mathfrak{t}}^{2}}, \ \mathcal{N}\mathfrak{u}\_{\mathfrak{t}} = \frac{\mathfrak{x}q\_{w}}{k\_{f}(T\_{w} - T\_{\infty})} \tag{15}$$

where the shear stress at wall *τw* and the constant surface heat flux *qw* are expressed as:

$$\pi\_{\overline{w}} = \mu\_{\imath f} \left( \frac{\partial u}{\partial y} \right)\_{y=0}, \ q\_{\overline{w}} = -k\_{\imath f} \left( \frac{\partial T}{\partial y} \right)\_{y=0} + (q\_r)\_{y=0}. \tag{16}$$

Substituting (9) into (16) and using (15), the following is obtained:

$$\operatorname{Re}\_x^{1/2} \mathbb{C}\_f = \frac{1}{(1-\phi)^{2.5}} f''(0),\\\operatorname{Re}\_x^{-1/2} N u\_\mathbb{X} = -\left(\frac{k\_{\rm nf}}{k\_f} + \frac{4}{3} Nr\right) \theta'(0). \tag{17}$$

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

The numerical results of the nonlinear ordinary differential equations given in Equations (10) and (11) together with the boundary conditions in Equation (12) indicates that for a particular range of the mixed convection parameter *λ*, there exist dual solutions (upper and lower branch solutions) for the various values of the selected governing parameters. Therefore, to validate which solution is in the stable flow, the stability of the dual solutions is tested by accommodating the stability analysis shown in Merkin [53]. To perform this, an unsteady form of the problem was considered. Equation (1) was retained, while Equations (2) and (7) were substituted by the following:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = u\_\epsilon \frac{du\_\epsilon}{dx} + \frac{\mu\_{nf}}{\rho\_{nf}} \frac{\partial^2 u}{\partial y^2} + \frac{(\rho \beta)\_{nf}}{\rho\_{nf}} (T - T\_\infty) \mathbf{g},\tag{18}$$

$$\frac{\partial u}{\partial t} + u \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial y} = a\_{nf} \frac{\partial^2 T}{\partial y^2} + \frac{16 \sigma^\* T\_{\infty}^3}{3 \left(\rho C\_p \right)\_{nf} k^\*} \frac{\partial^2 T}{\partial y^2} + \frac{Q\_0}{\left(\rho C\_p \right)\_{nf}} (T - T\_{\infty}), \tag{19}$$

where *t* represents the time. Analogous to the similarity transformation (9), the following new dimensionless functions *u*, *v* and *θ* have been introduced in conjunction to the similarity variable *η* which is the same as defined in (9), and the new similarity variable *τ* as follows:

$$u = a \text{tr} \frac{\partial f}{\partial \eta}, \ v = -\sqrt{a \nu\_f} f(\eta, \tau), \ \theta(\eta, \tau) = \frac{T - T\_{\infty}}{T\_w - T\_{\infty}}, \ \eta = \sqrt{\frac{a}{\nu\_f}} y, \ \tau = at. \tag{20}$$

Of note, with variables *u* and *v* given in the above, the equation of continuity (1) is identically satisfied. Next, after substituting the new similarity transformation (20) into Equations (18) and (19), we obtained the following equations:

$$\frac{1}{(1-\phi)^{2.5}}\frac{\partial^3 f}{\partial \eta^3} + \left(1-\phi+\phi\frac{\rho\_s}{\rho\_f}\right)\left(f\frac{\partial^2 f}{\partial \eta^2} + 1 - \left(\frac{\partial f}{\partial \eta}\right)^2 - \frac{\partial^2 f}{\partial \eta \partial \tau}\right) + \left(1-\phi+\phi\frac{(\rho\beta)\_s}{(\rho\beta)\_f}\right)\lambda\theta = 0,\tag{21}$$

$$\frac{1}{\Pr} \left( \frac{k\_{nf}}{k\_f} + \frac{4}{3} \text{Nr} \right) \frac{\partial^2 \theta}{\partial \eta^2} + \left( 1 - \phi + \phi \frac{\left( \rho \mathbb{C}\_p \right)\_s}{\left( \rho \mathbb{C}\_p \right)\_f} \right) \left( f \frac{\partial \theta}{\partial \eta} - \frac{\partial f}{\partial \eta} \theta + K\theta - \frac{\partial \theta}{\partial \tau} \right) = 0,\tag{22}$$

and were subjected to the boundary conditions:

$$f(\eta, \tau) = s, \,\frac{\partial f}{\partial \eta} = c, \theta(\eta, \tau) = 1 \text{ at } \eta = 0, \frac{\partial f}{\partial \eta} = 1, \,\theta(\eta, \tau) = 0 \text{ as } \eta \to \infty. \tag{23}$$

Next, to study the stability of the dual solutions, small disturbances of the growth (or decay) rate *γ* or better known as the unknown eigenvalue parameter, are taken in the form (see Weidman et al. [54]):

$$f(\eta, \tau) = f\_0(\eta) + \varepsilon^{-\gamma \tau} F(\eta, \tau), \; \theta(\eta, \tau) = \theta\_0(\eta) + \varepsilon^{-\gamma \tau} G(\eta, \tau), \tag{24}$$

where *f*0(*η*) and *<sup>θ</sup>*0(*η*) satisfied the problem (10)–(12). Besides, *<sup>F</sup>*(*η*, *<sup>τ</sup>*), *<sup>G</sup>*(*η*, *τ*) and all of the respective derivatives were assumed to be smaller when compared to *f*0(*η*), *<sup>θ</sup>*0(*η*) and its derivatives. By means of using (24), hence Equations (21) and (22) can be given as:

$$\frac{1}{(1-\phi)^{\frac{\alpha}{2}\delta}}\frac{\partial^{\frac{\alpha}{2}}F}{\partial\eta^{\frac{\alpha}{2}}}+\left(1-\phi+\phi\frac{\rho\_{s}}{\rho\_{f}}\right)\left(f\_{0}\frac{\partial^{2}F}{\partial\eta^{2}}+f\_{0}''F-2f\_{0}'\frac{\partial F}{\partial\eta}+\gamma\frac{\partial F}{\partial\eta}-\frac{\partial^{2}F}{\partial\eta\partial\tau}\right)+\left(1-\phi+\phi\frac{(\rho\delta)\_{x}}{(\rho\delta)\_{f}}\right)\lambda\theta = 0,\tag{25}$$

$$\rho\_{\rm TF}^{\frac{1}{2}} \left( \frac{k\_{\rm ref}}{k\_f} + \frac{4}{3} Nr \right) \frac{\partial^2 G}{\partial \eta^2} + \left( 1 - \phi + \phi \frac{\left( \rho C\_p \right)\_s}{\left( \rho C\_p \right)\_f} \right) \left( f\_0 \frac{\partial G}{\partial \eta} + F \theta\_0' - f\_0' G - \theta\_0 \frac{\partial F}{\partial \eta} + K G + \gamma G - \frac{\partial G}{\partial \overline{\tau}} \right) = 0,\tag{26}$$

together with the following boundary conditions:

$$\begin{array}{c} F(\eta, \tau) = 0, \ \frac{\partial F}{\partial \eta} = 0, G(\eta, \tau) = 0 \text{ at } \eta = 0, \\\ \frac{\partial F}{\partial \eta} = 0, G(\eta, \tau) = 0 \text{ as } \eta \to \infty. \end{array} \tag{27}$$

As proposed by Weidman et al. [54], the initial growth or decay of the solutions (24) is identified, by setting *τ* = 0, thus, giving *F* = *<sup>F</sup>*0(*η*) and *G* = *<sup>G</sup>*0(*η*). In this respect, the following linear eigenvalue problem was solved:

$$\frac{1}{(1-\phi)^{2.5}} F'''\_0 + \left(1 - \phi + \phi \frac{\rho\_s}{\rho\_f}\right) \left(f\_0 F''\_0 + f''\_0 F\_0 - 2f'\_0 F'\_0 + \gamma F'\_0\right) + \left(1 - \phi + \phi \frac{(\rho\beta)\_s}{(\rho\beta)\_f}\right) \lambda G\_0 = 0,\tag{28}$$

$$\frac{1}{\Pr} \left( \frac{k\_{nf}}{k\_f} + \frac{4}{3} Nr \right) \mathbf{G}\_0'' + \left( 1 - \phi + \phi \frac{\left( \rho \mathbf{C}\_p \right)\_s}{\left( \rho \mathbf{C}\_p \right)\_f} \right) \left( f\_0 \mathbf{G}\_0' + F\_0 \theta\_0' - f\_0' \mathbf{G}\_0 - \theta\_0 F\_0' + K \mathbf{G}\_0 + \gamma \mathbf{G}\_0 \right) = 0,\tag{29}$$

with the boundary conditions given by:

$$\begin{array}{l} F\_0(\eta) = 0, \; F\_0'(\eta) = 0, G\_0(\eta) = 0 \text{ at } \eta = 0, \\\ F\_0'(\eta) = 0, G\_0(\eta) = 0 \text{ as } \eta \to \infty. \end{array} \tag{30}$$

Indeed, it should be stated at this point, that the solutions *f*0(*η*) and *<sup>θ</sup>*0(*η*) were determined from the problem depicted in Equations (10)–(12). Upon obtaining the results, *f*0(*η*) and *<sup>θ</sup>*0(*η*) were again applied to Equations (28) and (29), and the linear eigenvalue problem (28)–(30) were solved. Harris et al. [55] proposed to relax a suitable boundary condition on *<sup>F</sup>*0(∞) = 0 or *<sup>G</sup>*0(∞) = 0 to determine a better range of *γ*. In the current study, the condition *<sup>F</sup>*0(∞) = 0 is relaxed and for a fixed value of *γ*, the linear eigenvalue problem (28)–(30) are solved, together with the new boundary condition; *F* 0 (0) = 1. Notably, it is worth mentioning that the solutions of the linear eigenvalue problem (28)–(30) provides an infinite set of eigenvalues *γ*1 < *γ*2 < *γ*3 < ... , where *γ*1 refers to the smallest eigenvalue. Furthermore, a positive *γ*1 reflects to an initial decay of disturbances and a stable flow. In contrast, a negative *γ*1 indicates an initial growth of disturbances and unstable flow.

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

The derived nonlinear ordinary differential equations given in Equations (10) and (11) along with the boundary conditions given in (12) were solved numerically and were obtained using the bvp4c programme in MATLAB (Matlab R2015a, MathWorks, Natick, MA, USA) for the selected values of the mixed convection parameter *λ*, solid volume fraction of nanoparticles *φ*, thermal radiation parameter *Nr*, heat source/sink parameter *K*, suction parameter *s* and stretching/shrinking parameter *c*. The range of *φ* values was taken as 0 ≤ *φ* ≤ 0.2, where *φ* = 0 indicates a regular base fluid, while the value of the Prandtl number was considered as Pr = 6.2 (water), except for comparisons with the prior case. The correlative output of the results obtained for *f* (0), with the ones obtained in Bachok et al. [56] for some values of *c* and *φ* with *λ* = *Nr* = *K* = *s* = 0 for Cu-water nanofluid, are presented in Table 2. Also, it was achieved that the present results were in very good alliance, which confirms that the numerical approach applied in this study is perfect, and therefore, the obtained results were believed to be accurate and correct.

The variations of the reduced skin friction coefficient *f* (0) and the reduced local Nusselt number −*θ*(0) against *λ* are shown in Figures 2–13 for several values of *φ*, *Nr*, *K*, *s* and *c*. It was observed that dual solutions (upper and lower branch solutions) occurred for Equations (10) and (11) subject to the boundary conditions (12) in the range of *λ* > *λ<sup>c</sup>*, where *λc* denotes the critical value of *λ*. Note that no solutions exist for *λ* < *λc* while a unique solution exists when *λ* = *λ<sup>c</sup>*. Also, it is obvious from these figures that the values of |*<sup>λ</sup>c*| increase as the parameters *φ*, *s* and *c* increase, therefore suggesting that these parameters widen the range of occurrence of dual solutions. Accordingly, it is further confirmed that the presence of nanoparticles, heat sink, suction and stretching sheet could decelerate the separation of the boundary layer, while the presence of thermal radiation, heat source and shrinking sheet could accelerate the separation of the boundary layer. Also, the values of −*θ*(0) are always positive for the upper branch solution which is due to the heat being transferred from the

hot surface of the stretching/shrinking sheet to the cold fluid. The reverse trend is observed in the case of the lower branch solution, i.e., −*θ* (0) becomes unbounded as *λ* → 0<sup>+</sup> and *λ* → 0− .


**Table 2.** Values of *f* (0) for various values of the stretching/shrinking parameter *c* and the solid volume fraction of nanoparticles *φ* for Cu-water nanofluid.

"[ ]" refers to the lower branch solution.

**Figure 2.** Variation of the reduced skin friction coefficient *f* (0) with the mixed convection parameter *λ* for several values of the volume fraction of nanoparticles *φ* when *Nr* = 0.1, *K* = −0.1, *c* = −1 and *s* = 1 for Cu-water nanofluid.

**Figure 3.** Variation of the reduced local Nusselt number −*θ* (0) with the mixed convection parameter *λ* for several values of the volume fraction of nanoparticles *φ* when *Nr* = 0.1, *K* = −0.1, *c* = −1 and *s* = 1 for Cu-water nanofluid.

**Figure 4.** Variation of the reduced skin friction coefficient *f* (0) with the mixed convection parameter *λ* for several values of the radiation parameter *Nr* when *φ* = 0.01, *K* = −0.1, *s* = 1 and *c* = −1 for Cu-water nanofluid.

**Figure 5.** Variation of the reduced local Nusselt number −*θ* (0) with the mixed convection parameter *λ* for several values of the radiation parameter *Nr* when *φ* = 0.01, *K* = −0.1, *s* = 1 and *c* = −1 for Cu-water nanofluid.

**Figure 6.** Variation of the reduced skin friction coefficient *f* (0) with the mixed convection parameter *λ* for various values of the heat source/sink parameter *K* when *φ* = 0.01, *Nr* = 0.1, *s* = 1 and *c* = −1 for Cu-water nanofluid.

**Figure 7.** Variation of the reduced local Nusselt number −*θ*(0) with the mixed convection parameter *λ* for various values of the heat source/sink parameter *K* when *φ* = 0.01, *Nr* = 0.1, *s* = 1 and *c* = −1 for Cu-water nanofluid.

**Figure 8.** Variation of the reduced skin friction coefficient of *f* (0) with mixed convection parameter *λ* for various values of the suction parameter *s* when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *c* = −1 for Cu-water nanofluid.

**Figure 9.** Variation of the reduced local Nusselt number −*θ*(0) with the mixed convection parameter *λ* for various values of the suction parameter *s* when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *c* = −1 for Cu-water nanofluid.

**Figure 10.** Variation of the reduced skin friction coefficient of *f* (0) with the mixed convection parameter *λ* for several values of *c* (shrinking sheet) when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *s* = 1 for Cu-water nanofluid.

**Figure 11.** Variation of the reduced local Nusselt number −*θ*(0) with the mixed convection parameter *λ* for several values of *c* (shrinking sheet) when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *s* = 1 for Cu-water nanofluid.

**Figure 12.** Variation of the reduced skin friction coefficient of *f* (0) with the mixed convection parameter *λ* for several values of *c* (stretching sheet) when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *s* = 1 for Cu-water nanofluid.

**Figure 13.** Variation of the reduced local Nusselt number −*θ*(0) with the mixed convection parameter *λ* for various values of *c* (stretching sheet) when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1 and *s* = 1 for Cu-water nanofluid.

The influences of the solid volume fraction of nanoparticles *φ* on *f* (0) and −*θ*(0) are illustrated in Figures 2 and 3. Figure 2 illustrates that *f* (0) increases smoothly with an increasing value of *φ*. Brinkman [57] explained that an increment in the nanoparticle volume fraction increases the fluid's viscosity. Hence, this situation contributed towards increasing the skin friction along the surface. Figure 3 depicts that as the volume fraction of nanoparticles increases which consequently decreases the rate of heat transfer at the surface of <sup>−</sup>*θ*(0). This occurs because the nanoparticles increase the viscosity, density and conductivity. However, they may decrease the heat capacitance of the nanofluid *<sup>ρ</sup>Cpn f* (see MacDevette et al. [58]). Therefore, there is a trade-off between the enhanced properties, increased viscosity, decreased *<sup>ρ</sup>Cpnf*and possible decrease in the heat transfer coefficient.

Figures 4 and 5 display the impacts of the thermal radiation parameter *Nr* on *f* (0) and <sup>−</sup>*θ*(0), respectively. Figure 4 also illustrates that the reduced skin friction coefficient increases as *Nr* is reduced in the case of the opposing flow. The reverse trend is noted in the event of the assisting flow and an increment in thermal radiation enhances the transmission of energy between the nanoparticles on a heated surface (assisting flow). Eventually, the thickness of the momentum boundary layer becomes thinner and increases the wall shear stress which raise the values of *f* (0). The decrement trend in the rate of heat transfer at the surface leads to an increment in the thermal radiation. Notably, this is in accordance with the result shown in Figure 5 where the reduced local Nusselt number −*θ*(0) decreases with increasing *Nr*. Notwithstanding, this fact can be explained as follows. As the influence of the thermal radiation becomes stronger, the thermal boundary layer thickness increases and further decreases the values of the rate of heat transfer. Therefore, the usage of nanofluids having thermal radiation cannot improve the cooling of the heated sheet.

Figure 6 displays the variations of *f* (0) for various values of the heat source/sink parameter *K*. It is noticeable that the decrement of *f* (0) as the heat source/sink parameter increases from the negative (heat sink) to the positive values (heat source) in the event of the opposing flow. A contrary trend is noticeable in the event of the assisting flow. The influence of the heat source/sink parameter *K* on the reduced local Nusselt number −*θ*(0) is illustrated in Figure 7. The results presented in this figure indicate that the rate of heat transfer at the surface decreases as the heat source/sink parameter increases from negative (heat sink) to positive values (heat source). Indeed, this is because the higher heat source effect can increase the thermal boundary layer thickness that reduces the rate of heat transfer.

Figure 8 illustrates the effect of the suction on the reduced skin friction coefficient. Further, it is observed that the values of *f* (0) increase with the suction. This is because the influence of suction at the boundary that slows down the nanofluid motion and increases the velocity gradient at the surface. By observing Figure 9, it is evident that the values of −*θ*(0) that represents the heat transfer rate at the

surface, increases with the suction. Precisely, the increase in the magnitude of the suction parameter consequently increased the rate of heat transfer. This is because the increasing suction decreases the thermal boundary layer thickness and in return increases the temperature gradient at the surface.

Figures 10–13 are shown to present the impacts of stretching/shrinking parameter *c* on *f* (0) and −*θ* (0). Notably, it is evident from these figures that *f* (0) and −*θ* (0) were higher in the event of the stretching sheet, in comparison with the case of the shrinking sheet. Hence, the stretching parameter provided the most significant effects upon the skin friction along the surface and the heat flux at the surface. Indeed, this suggests that the stretching sheet enhances the rate of heat transfer, whereas the shrinking sheet inhibits the effect of the heat transfer rate. Also, the increment in the value of the stretching parameter further increases the effect of free convection.

The variations of *f* (0) and −*θ* (0) for three distinct types of nanofluids with nanoparticles containing Cu, Al2O3 and TiO2 are shown in Figures 14 and 15 From the figures, it is evident that there is little difference in the value of the reduced skin friction coefficient and the reduced local Nusselt number for TiO2- and Al2O3-water nanofluids. Further, it is also discovered that the values of *f* (0) and −*θ* (0) are highest for Cu-water nanofluid, followed by TiO2- and Al2O3-water nanofluids which is due to Cu having the highest value of thermal conductivity in comparison to the other nanoparticles.

**Figure 14.** Variation of the reduced skin friction coefficient of *f* (0) with the mixed convection parameter *λ* when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1, *s* = 1 and *c* = −1 for Cu-water, Al2O3-water and TiO2-water nanofluids.

**Figure 15.** Variation of the reduced local Nusselt number −*θ* (0) with the mixed convection parameter *λ* when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1, *s* = 1 and *c* = −1 for Cu-water, Al2O3-water and TiO2-water nanofluids.

The velocity and temperature profiles for Cu-water nanofluid with different values of the volume fraction of the nanoparticles *φ* are displayed in Figures 16 and 17. It was discovered that an increment in the value of *φ*, increased the velocity of the fluid. Since an increment in *φ* is believed to increase the fluid's viscosity as suggested by Brinkman [57], hence it increases the skin friction along the permeable vertical shrinking flat plate. This statement can be proved by the velocity profiles as shown in Figure 16, as an increment in *φ* reduces the momentum boundary layer thickness. Furthermore, from Figure 17, we can see that the temperature of the fluid increases with an increase in *φ*, therefore, suggesting that the temperature of the nanofluids can be managed by increasing or decreasing the volume fraction of the nanoparticles in the base fluid.

**Figure 16.** Velocity profiles *f* (*η*) for various values of the volume fraction of the nanoparticles *φ* when *Nr* = 0.1, *K* = 0.1, *s* = 1, *c* = −1 and *λ* = −2 (opposing flow) for Cu-water nanofluid.

**Figure 17.** Temperature profiles *<sup>θ</sup>*(*η*) for various values of the volume fraction of the nanoparticles *φ* when *Nr* = 0.1, *K* = 0.1, *s* = 1, *c* = −1 and *λ* = −2 (opposing flow) for Cu-water nanofluid.

The velocity and temperature profiles for Cu-, Al2O3- and TiO2-water nanofluids are shown in Figures 18 and 19, which indicate that by utilizing distinct types of nanofluids, the values of velocity and temperature change. Furthermore, we discovered that Cu-water nanofluid has higher velocity distribution and lower temperature distribution in comparison to the other two nanofluids for the upper branch solution. Also, it is discovered that the velocity and temperature profiles for Al2O3-water and TiO2-water nanofluids nearly coincide with each other. Besides, the Cu-water nanofluid (compared

to Al2O3-water and TiO2-water nanofluids) has thinner momentum and thermal boundary layer thickness which is due to the fact Cu nanoparticles have the highest thermal conductivity value in comparison to the other two kinds of nanoparticles. Accordingly, the reduced value of thermal diffusivity causes higher temperature gradients and therefore, the higher enhancement in heat transfers. Notwithstanding, the Cu nanoparticle possess high values of thermal diffusivity, and therefore, lowers the temperature gradients that affects the performance of Cu nanoparticle.

**Figure 18.** Velocity profiles *f* (*η*) for Cu-water, Al2O3-water and TiO2-water nanofluids when *φ* = 0.01, *Nr* = 0.1, *K* = 0.1, *s* = 1, *c* = −1 and *λ* = −2 (opposing flow).

**Figure 19.** Temperature profiles *<sup>θ</sup>*(*η*) for Cu-water, Al2O3-water and TiO2-water nanofluids when *φ*= 0.01, *Nr* = 0.1, *K* = 0.1, *s* = 1, *c* = −1 and *λ* = −2 (opposing flow).

From observing Figures 16–19, is apparent that the velocity and temperature profiles for both the upper and lower branch solutions satisfied the far field boundary conditions (12) asymptotically. Thus, aids in validating and supporting the numerical results retrieved for the boundary value problem (10)–(12) and proved the occurrence of the dual nature of the solutions as shown in Figures 2–15. Also, the lower branch solution for the corresponding profiles displayed a larger boundary layer thickness when compared to the upper branch solution.

A stability analysis was undertaken to test the stability of the dual solutions. Hence, the smallest eigenvalue *γ*1 was determined by solving the linear eigenvalue problem (28)–(30) using the bvp4c programme in MATLAB. Table 3 depicts the smallest eigenvalues for Cu-water nanofluid for some

values of parameters *K*, *Nr* and *λ* when *s* = 1 and *c* = −1. The results in Table 3 further indicate that *γ*1 is negative for the lower branch solution, while *γ*1 is positive for the upper branch solution due to their correspondence to the initial decay of disturbances, thus signifying a stable flow. Furthermore, for both the upper and lower branch solutions, the values of |*<sup>γ</sup>*1| → 0 as *λ* approaches *λ<sup>c</sup>*, is consistent with Merkin [53], Weidman et al. [54] and Harris et al. [55].


**Table 3.** Smallest eigenvalues *γ*1 for various values of *K*, *Nr* and *λ*.

## **5. Conclusions**

In this paper, the mixed convection flow near a stagnation-point over a vertical stretching/shrinking sheet with suction, thermal radiation and heat source/sink in Cu-, Al2O3- and TiO2-water nanofluids has been investigated numerically. The main limitation of this study is the mathematical nanofluid model used which considers boundary layer approximations, thus, we are unable to ge<sup>t</sup> the dual solutions beyond the separation point of a laminar boundary layer by utilizing the boundary layer approximations. To achieve the dual solutions beyond this point, the full Navier-Stokes equations are solved. Unfortunately, this case is way beyond the scope of the present research. Another limitation is that we utilized the mathematical nanofluid model suggested by Tiwari and Das [49], which does not consider the impacts of the Brownian motion and thermophoresis on nanofluids. It is worth acknowledging that a mathematical nanofluid model that refers to the Brownian motion and thermophoresis has been suggested by Buongiorno [59]. The results from this study may be summarised as follows:


It is observed from Table 2 that the values of *f* (0) for *φ* = 0 (classical viscous fluid) are lower than those for a nanofluid (*φ* = <sup>0</sup>). *f* (0) increases almost monotonically with increasing the nanoparticle volume fraction *φ*.

**Author Contributions:** Numerical analysis were performed by A.J. and R.N. These authors also explained the results and wrote the manuscript. The literature review was written by I.P. and he co-wrote the manuscript. All authors originated and developed the problem and reviewed the manuscript.

**Funding:** This research was funded by the research university gran<sup>t</sup> from the Universiti Kebangsaan Malaysia (UKM), gran<sup>t</sup> number DIP-2017-009. The work of Ioan Pop was supported from Unitatea Executivă pentru Finan¸tarea Învă¸tământului Superior, a Cercetării, Dezvoltării ¸si Inovării (UEFISCDI), Romania, gran<sup>t</sup> number PN-III-P4-ID-PCE-2016-0036.

**Acknowledgments:** The authors would like to thank the reviewers for their constructive comments which beautify the quality of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
