2.1. Model

Consider the flow of a nanofluid over a semi-infinite flat plate with free stream velocity $U_{\infty }^*$ . (Here, an asterisk denotes dimensional quantities.) The model is given in Cartesian coordinates $\boldsymbol{x}^*=(x^*,y^*,z^*)$ , where $x^*$ measures the distance along the surface of the flat plate, $y^*$ denotes the direction normal to the plate and $z^*$ the spanwise direction. Consequently, the governing system of equations comprise the continuity, momentum and energy equations for fluid motion (Ruban & Gajjar Reference Ruban and Gajjar2014), along with a continuity equation for the nanoparticles (Buongiorno Reference Buongiorno2006; Avramenko et al. Reference Avramenko, Blinov and Shevchuk2011; MacDevette et al. Reference MacDevette, Myers and Wetton2014),

(2.1a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \frac {\partial \rho ^*}{\partial t^*} + {\nabla} ^*\boldsymbol{\cdot }(\rho ^*\boldsymbol{u}^*) = 0, \end{align}

(2.1b) \begin{align}& \rho ^*\left (\frac {\partial \boldsymbol{u}^*}{\partial t^*} + (\boldsymbol{u}^*\boldsymbol{\cdot }{\nabla} ^*)\boldsymbol{u}^*\right ) = -{\nabla} ^*p^* + {\nabla} ^*\boldsymbol{\cdot }\left (\mu ^*\left ({\nabla} ^*\boldsymbol{u}^* + \left ({\nabla} ^*\boldsymbol{u}^*\right )^T - \frac {2}{3}{\nabla} ^*\boldsymbol{\cdot }\boldsymbol{u}^*\unicode{x1D644}\right )\right )\!, \end{align}

(2.1c) \begin{align}& \rho ^*\left (\frac {\partial \left (c^*T^*\right )}{\partial t^*} + (\boldsymbol{u}^*\boldsymbol{\cdot }{\nabla} ^*)\left (c^*T^*\right )\right ) = {\nabla} ^*\boldsymbol{\cdot }\left (k^*{\nabla} ^* T^*\right ) \nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad + (\rho ^*c^*)_{\textit{np}}\left (D_B^*{\nabla} ^*\phi +D_T^*\frac {{\nabla} ^*T^*}{T^*}\right )\boldsymbol{\cdot }{\nabla} ^*T^*, \end{align}

(2.1d) \begin{align}&\qquad\qquad\qquad \frac {\partial \phi }{\partial t^*}+{\nabla} ^*\boldsymbol{\cdot }(\phi \boldsymbol{u}^*) = {\nabla} ^*\boldsymbol{\cdot }\left (D_B^*{\nabla} ^*\phi + D_T^*\frac {{\nabla} ^*T^*}{T^*}\right )\!, \end{align}

for a velocity $\boldsymbol{u}^*=(u^*,v^*,w^*)$ , pressure $p^*$ , temperature $T^*$ and dimensionless nanoparticle volume concentration $\phi$ . Here, $\unicode{x1D644}$ is the identity matrix.

The density of the nanofluid $\rho ^*$ is defined using the law of mixtures as

(2.2) \begin{align} \rho ^* = (1-\phi )\rho _{\textit{bf}}^* + \phi \rho _{\textit{np}}^*, \end{align}

where subscripts $\textit{bf}$ and $np$ represent quantities associated with the base fluid and nanoparticles, respectively. In addition,

(2.3) \begin{align} \rho ^*c^* = (1-\phi )(\rho ^*c^*)_{\textit{bf}} + \phi (\rho ^*c^*)_{\textit{np}}, \end{align}

where $c^*$ denotes the specific heat capacity of the nanofluid, while the thermal conductivity of the nanofluid $k^*$ is given by the Maxwell (Reference Maxwell1881) model:

(2.4) \begin{align} k^* = \left (\frac {k_{\textit{np}}^*+2k_{\textit{bf}}^*+2\phi \left(k_{\textit{np}}^*-k_{\textit{bf}}^*\right)}{k_{\textit{np}}^*+2k_{\textit{bf}}^*-\phi \left(k_{\textit{np}}^*-k_{\textit{bf}}^*\right)}\right )k_{\textit{bf}}^*. \end{align}

Alternative models for $k^*$ may be considered as discussed in Wang & Mujumdar (Reference Wang and Mujumdar2008a ).

The dynamic viscosity of the nanofluid $\mu ^*$ , used throughout the subsequent study, is given by the Brinkman (Reference Brinkman1952) model,

(2.5) \begin{align} \mu ^* = \frac {\mu _{\textit{bf}}^*}{(1-\phi )^{2.5}}, \end{align}

for a base fluid dynamic viscosity $\mu _{\textit{bf}}^*$ . The Brinkman relation is known to under predict the dynamic viscosity for $\phi \gt 0.01$ (MacDevette et al. Reference MacDevette, Myers and Wetton2014). However, for theoretical purposes and to demonstrate trends, here we consider nanoparticle volume concentrations $\phi$ up to 10 % of the fluid volume. Similar to the thermal conductivity $k^*$ , alternative models may be considered for the dynamic viscosity $\mu ^*$ , as listed in Wang & Mujumdar (Reference Wang and Mujumdar2008a ) and Mishra et al. (Reference Mishra, Mukherjee, Nayak and Panda2014), which encompass properties such as the size and distribution of nanoparticles. For instance, Batchelor (Reference Batchelor1977) modelled the dynamic viscosity as

(2.6a) \begin{align} \mu ^* = \mu _{\textit{bf}}^*\left(1+2.5\phi +6.2\phi ^2\right)\!, \end{align}

whereas Pak & Cho (Reference Pak and Cho1998) and Maiga et al. (Reference Maiga, Nguyen, Galanis and Roy2004) obtained the correlations

(2.6b,c) \begin{align} \mu ^* = \mu _{\textit{bf}}^*\left(1+39.11\phi +533.9\phi ^2\right) \quad \textrm {and} \quad \mu ^* = \mu _{\textit{bf}}^*\left(1+7.3\phi +123\phi ^2\right)\!, \end{align}

for nanofluids inside circular pipes and tubes, respectively.

The latter two terms of (2.1c ) and the two terms on the right-hand side of (2.1d ) model the respective effects of BM and TP, with coefficients

(2.7a,b) \begin{align} D_B^* = \frac {k_B^*T^*}{3\pi \mu _{\textit{bf}}^*d_{\textit{np}}^*}\equiv C_B^*T^* \quad \textrm {and} \quad D_T^* = \frac {\beta _T\mu _{\textit{bf}}^*\phi }{\rho _{\textit{bf}}^*}\equiv C_T^*\phi . \end{align}

Here, $k_B^*$ denotes the Boltzmann constant, $d_{\textit{np}}^*$ the diameter of the nanoparticles and the proportionality constant

(2.8) \begin{align} \beta _T = 0.26\frac {k_{\textit{bf}}^*}{2k_{\textit{bf}}^*+k_{\textit{np}}^*}, \end{align}

as given in McNab & Meisen (Reference McNab and Meisen1973), Buongiorno (Reference Buongiorno2006) and MacDevette et al. (Reference MacDevette, Myers and Wetton2014).

The nanofluid flow is subject to the no-slip condition and the fixed temperature condition on the plate surface

(2.9a,b) \begin{align} \boldsymbol{u}^*=0 \quad \textrm {and} \quad T^*=T_w^* \quad \textrm {on} \quad y^*=0, \end{align}

where $T_w^*$ denotes the constant wall temperature. (Here, a subscript $w$ references wall conditions.) In addition,

(2.9c) \begin{align} D_B^*\frac {\partial \phi }{\partial y^*}+\frac {D_T^*}{T^*}{\frac {\partial T^*}{\partial y^*}}=0 \quad \textrm {on} \quad y^*=0, \end{align}

following Avramenko et al. (Reference Avramenko, Blinov and Shevchuk2011), which imposes that the total flux of nanoparticles at the plate surface is zero. Finally, in the far-field, the flow is subject to the free stream conditions

(2.10a–f) \begin{align} u^*&{}\rightarrow {} U_{\infty }^*, \qquad &v^*{}\rightarrow {}& 0, \qquad &w^*{}\rightarrow {}& 0, \nonumber \\ p^*&{}\rightarrow {} p_{\infty }^*, \qquad &T^*{}\rightarrow {}& T_{\infty }^*, \qquad &\phi {}\rightarrow {}&\phi _{\infty } \qquad \textrm {as} \quad y^*\rightarrow \infty , \end{align}

where $p_{\infty }^*$ , $T_{\infty }^*$ and $\phi _{\infty }$ denote the free stream pressure, the free stream temperature and the dimensionless free stream nanoparticle volume concentration, respectively. Figure 1 shows a schematic diagram of the nanofluid flow over a flat plate.

Figure 1.

Diagram of a nanofluid flow, composed of a base fluid ( $\textit{bf}$ ) and nanoparticles ( $np$ ) over a flat plate. Here, $\delta ^*$ represents the boundary-layer thickness.

2.2. Non-dimensionalisation

The governing system of equations (2.1) are non-dimensionalised by setting

(2.11a–i) \begin{align} \boldsymbol{x}^*&{}={}L^*\boldsymbol{x}, \qquad &\boldsymbol{u}^*{}={}&U_{\infty }^*\boldsymbol{u}, \qquad &t^*{}={}&L^*t/U_{\infty }^*, \nonumber \\ p^*&{}={}p_{\infty }^*+\rho _{\textit{bf}}^*U_{\infty }^{*2}p, \qquad &T^*{}={}&T_{\infty }^*T, \qquad &\rho ^*{}={}&\rho _{\textit{bf}}^*\rho , \nonumber \\ \mu ^*&{}={}\mu _{\textit{bf}}^*\mu , \qquad &c^*{}={}&c_{\textit{bf}}^*c, \qquad &k^*{}={}&k_{\textit{bf}}^*k, \end{align}

for a characteristic length scale $L^*$ . Consequently, (2.1) becomes

(2.12a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\rho \boldsymbol{u}\right ) = 0, \end{align}

(2.12b) \begin{align}&\qquad \rho \left (\frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}\right ) = -\boldsymbol{\nabla }p + \frac {1}{\textit{Re}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\mu \left (\boldsymbol{\nabla }\boldsymbol{u} + \left (\boldsymbol{\nabla }\boldsymbol{u}\right )^T - \frac {2}{3}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}\unicode{x1D644}\right )\right )\!, \end{align}

(2.12c) \begin{align}& \rho \left (\frac {\partial \left (cT\right )}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\left (cT\right )\right ) = \frac {1}{\textit{Re} \textit{Pr}}\boldsymbol{\nabla }\boldsymbol{\cdot }(k\boldsymbol{\nabla }T) + \frac {1}{\textit{Re} \textit{Pr} \textit{Le}} \left ( T\boldsymbol{\nabla }\phi + \frac {\phi \boldsymbol{\nabla }T}{N_{{\textit{BT}}}T} \right )\boldsymbol{\cdot }\boldsymbol{\nabla }T, \end{align}

(2.12d) \begin{align}&\qquad\qquad\qquad \frac {\partial \phi }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{u}) = \frac {1}{\textit{Re} \textit{Sc}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (T\boldsymbol{\nabla }\phi +\frac {\phi \boldsymbol{\nabla }T}{N_{ {\textit{BT}}}T}\right )\!, \end{align}

where

(2.13a) \begin{align} \rho &= 1+(\hat {\rho }-1)\phi \qquad \qquad \,\,\,\, \textrm {for} \quad \hat {\rho } = \frac {\rho _{\textit{np}}^*}{\rho _{\textit{bf}}^*}, \end{align}

(2.13b) \begin{align} \rho c &= 1+(\hat {\rho }\hat {c}-1)\phi \qquad \qquad \,\; \textrm {for} \quad \hat {\rho }\hat {c} = \frac {(\rho ^*c^*)_{\textit{np}}}{(\rho ^*c^*)_{\textit{bf}}}, \end{align}

(2.13c) \begin{align} k &= \left (\frac {\hat {k}+2+2(\hat {k}-1)\phi }{\hat {k}+2-(\hat {k}-1)\phi }\right ) \quad \textrm {for} \quad \hat {k} = \frac {k_{\textit{np}}^*}{k_{\textit{bf}}^*}. \end{align}

Moreover, in the case of the Brinkman (Reference Brinkman1952) viscosity model, given by (2.5), the non-dimensional dynamic viscosity is given as

(2.14) \begin{align} \mu = \frac {1}{(1-\phi )^{2.5}}. \end{align}

Similar representations for $\mu$ are given for the Batchelor (Reference Batchelor1977), Pak & Cho (Reference Pak and Cho1998) and Maiga et al. (Reference Maiga, Nguyen, Galanis and Roy2004) models.

Table 1.

Thermophysical properties of water and various materials used for nanoparticles, as given in Buongiorno (Reference Buongiorno2006), Wang & Mujumdar (Reference Wang and Mujumdar2008a ), Bachok et al. (Reference Bachok, Ishak and Pop2011), MacDevette et al. (Reference MacDevette, Myers and Wetton2014), Turkyilmazoglu (Reference Turkyilmazoglu2014, Reference Turkyilmazoglu2020) and at https://periodictable.com/Elements. Here, free stream temperature $T_{\infty }^*=300$ K, nanoparticle diameter $d_{\textit{np}}^*=20$ nm and Prandtl number ${\textit{Pr}}=6.85$ . The ratios $\hat {\rho }$ , $\hat {k}$ and $\hat {c}$ are based on water as the base fluid.

Figure 2 compares the four models of the non-dimensional dynamic viscosity $\mu$ along with the non-dimensional density $\rho$ , thermal conductivity $k$ and specific heat capacity $c$ for Cu nanoparticles in water (see table 1 for thermophysical properties). These quantities are plotted as functions of the free stream nanoparticle volume concentration $\phi _{\infty }$ . As $\phi _{\infty }$ increases, the Brinkman and Batchelor viscosity models show a similar rate of increase, while the Pak–Cho and Maiga viscosity models exhibit a more rapid increase. In addition, $\rho$ also increases with $\phi _{\infty }$ . Furthermore, $k$ increases, improving the flows heat transfer capability, while $c$ exhibits a reduction, causing temperature changes within the flow to occur more rapidly.

Figure 2.

( $a$ ) Non-dimensional dynamic viscosity $\mu$ as a function of $\phi _{\infty }$ , for the Brinkman (Reference Brinkman1952), Batchelor (Reference Batchelor1977), Pak & Cho (Reference Pak and Cho1998) and Maiga et al. (Reference Maiga, Nguyen, Galanis and Roy2004) models. ( $b$ ) Non-dimensional density $\rho$ , specific heat capacity $c$ and thermal conductivity $k$ as a function of $\phi _{\infty }$ , for Cu nanoparticles in water. Refer to table 1 for fluid and nanoparticle properties.

The dimensionless Reynolds, Prandtl, Lewis and Schmidt numbers are defined as

(2.15a,b) \begin{align} \textit{Re} &{}={} \frac {U_{\infty }^*L^*\rho _{\textit{bf}}^*}{\mu _{\textit{bf}}^*}, \quad &\textit{Pr} {}={}& \frac {\mu _{\textit{bf}}^* c_{\textit{bf}}^*}{k_{\textit{bf}}^*}, \end{align}

(2.15c,d) \begin{align} \textit{Le} &{}={} \frac {k_{\textit{bf}}^*}{(\rho ^*c^*)_{\textit{np}}C_B^*T_{\infty }^*}, \quad &{} \textit{Sc} ={}& \frac {\mu _{\textit{bf}}^*}{\rho _{\textit{bf}}^*C_B^*T_{\infty }^*}, \end{align}

while the ratio of BM to TP is given as

(2.16) \begin{align} \quad N_{{\textit{BT}}}=\frac {C_B^*T_{\infty }^*}{C_T^*}. \end{align}

Finally, the boundary conditions (2.9) on the plate surface are recast as

(2.17a,b) \begin{align} \boldsymbol{u}=0 \quad \textrm {and} \quad T = T_w \left (\equiv \frac {T_w^*}{T_{\infty }^*}\right ) \quad \textrm {on} \quad y=0, \end{align}

and

(2.17c) \begin{align} T\frac {\partial \phi }{\partial y}+\frac {\phi }{N_{{\textit{BT}}}T}{\frac {\partial T}{\partial y}}=0 \quad \textrm {on} \quad y=0, \end{align}

while the boundary conditions (2.10) in the free stream are given as

(2.18a–f) \begin{align} u&{}\rightarrow {} 1, \qquad &v{}\rightarrow {}& 0, \qquad &w{}\rightarrow {}& 0, \nonumber \\ p&{}\rightarrow {} 0, \qquad &T{}\rightarrow {}& 1, \qquad &\phi {}\rightarrow {}&\phi _{\infty } \qquad \textrm {as} \quad y\rightarrow \infty . \end{align}

Table 1 presents the thermophysical properties of various materials used for nanoparticles. The non-dimensional ratios $\hat {\rho }$ , $\hat {k}$ and $\hat {c}$ are based on water as the base fluid, where the Prandtl number ${\textit{Pr}}=6.85$ , while the Lewis number $Le$ , the Schmidt number ${\textit{Sc}}$ and the ratio $N_{{\textit{BT}}}$ are given for the free stream temperature $T_{\infty }^*=300$ K and the nanoparticle diameter $d_{\textit{np}}^*=20$ nm. Both $Le$ and ${\textit{Sc}}$ are of the order $10^4$ for all materials listed in table 1.

3.1. Boundary-layer equations

Following the derivation of Ruban (Reference Ruban2017), the steady, two-dimensional boundary-layer equations are obtained by assuming a zero pressure gradient, setting $w=0$ , and considering solutions that are independent of the $z$ -direction and time $t$ . On introducing the Prandtl boundary-layer transformation

(3.1) \begin{align} y = {\textit{Re}}^{-1/2}Y, \end{align}

with

(3.2a–h) \begin{align} u(x,y) &= U_{\!B}(x,Y), \quad v(x,y) {}={} {\textit{Re}}^{-1/2}V_{\!B}(x,Y), \nonumber \\ T(x,y) &= T_{\!B}(x,Y), \quad \phi (x,y) {}={} \phi _B(x,Y), \nonumber \\ \mu (x,y) &= \mu _B(x,Y), \quad\! \rho (x,y) {}={} \rho _B(x,Y), \nonumber \\ c(x,y) &= c_{\!B}(x,Y), \quad\, k(x,y) {}={} k_B(x,Y), \end{align}

and letting ${\textit{Re}}\rightarrow \infty$ , the non-dimensional governing equations (2.12) become

(3.3a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \frac {\partial (\rho _B U_{\!B})}{\partial x} + \frac {\partial (\rho _B V_{\!B})}{\partial Y} = 0, \end{align}

(3.3b) \begin{align}&\qquad\qquad\qquad\qquad \rho _B\left (U_{\!B}\frac {\partial U_{\!B}}{\partial x} + V_{\!B}\frac {\partial U_{\!B}}{\partial Y}\right ) = \frac {\partial }{\partial Y}\left (\mu _B\frac {\partial U_{\!B}}{\partial Y}\right )\!, \end{align}

(3.3c) \begin{align}&\quad \rho _B \left (U_{\!B}\frac {\partial (c_{\!B}T_{\!B})}{\partial x}+V_{\!B}\frac {\partial (c_{\!B}T_{\!B})}{\partial Y}\right ) = \frac {1}{\textit{Pr}}\frac {\partial }{\partial Y}\left (k_B\frac {\partial T_{\!B}}{\partial Y}\right ) \nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad + \frac {1}{\textit{Pr Le}} \left ( T_{\!B}\frac {\partial \phi _B}{\partial Y}\frac {\partial T_{\!B}}{\partial Y} +\frac {\phi _B}{N_{{\textit{BT}}}T_{\!B}}\left (\frac {\partial T_{\!B}}{\partial Y}\right )^2\right )\!, \end{align}

(3.3d) \begin{align}&\qquad\qquad\qquad \frac {\partial (\phi _B U_{\!B})}{\partial x} + \frac {\partial (\phi _B V_{\!B})}{\partial Y} = \frac {1}{\textit{Sc}}\frac {\partial }{\partial Y}\left (T_{\!B}\frac {\partial \phi _B}{\partial Y}+\frac {\phi _B}{N_{{\textit{BT}}}T_{\!B}}\frac {\partial T_{\!B}}{\partial Y}\right ). \end{align}

A self-similar solution is then sought using the similarity variable $\eta =Y/\sqrt {x}$ , coupled with the Dorodnitsyn–Howarth transformation

(3.4) \begin{align} \xi =\int _0^{\eta }\rho (\grave {\eta }) \; \textrm {d}\grave {\eta }, \end{align}

with

(3.5a–h) \begin{align} U_{\!B}(x,Y)&{}={}f'(\xi ), \quad V_{\!B}(x,Y){}={}\frac {1}{2\sqrt {x}}\left (\eta f'-\frac {f}{\rho }\right )\!, \nonumber \\ T_{\!B}(x,Y)&{}={}\theta (\xi ), \quad\,\, \phi _B(x,Y){}={}\varphi (\xi ), \nonumber \\ \mu _B(x,Y)&{}={}\mu (\xi ), \quad\, \rho _B(x,Y){}={}\rho (\xi ), \nonumber \\ c_{\!B}(x,Y)&{}={}c(\xi ), \quad\,\, k_B(x,Y){}={}k(\xi ). \end{align}

(For notational simplicity, $\mu$ , $\rho$ , $c$ and $k$ are reused to denote their similarity profiles.) Consequently, the following boundary-layer equations are derived:

(3.6a) \begin{align} 2\left (\rho \mu f''\right )'+ff''=0, \end{align}

(3.6b) \begin{align} 2\left (\rho k \theta '\right )' + Pr f\left (c\theta \right )' + \frac {2\rho \theta '}{Le} \left ( \theta \varphi ' + \frac {\varphi \theta '}{N_{{\textit{BT}}}\theta }\right ) = 0, \end{align}

(3.6c) \begin{align} \frac {2\rho ^2}{\textit{Sc}}\left (\rho \left (\theta \varphi ' + \frac {\varphi \theta '}{N_{{\textit{BT}}}\theta }\right )\right )^{\!\prime}+f\varphi '=0, \end{align}

subject to the boundary conditions

(3.6d–f) \begin{align} f=f'=0, \; \theta =T_w \quad \textrm {on} \quad \xi &=0, \end{align}

(3.6g) \begin{align} \theta \varphi '+\frac {\varphi \theta '}{N_{{\textit{BT}}}\theta }=0 \quad \textrm {on} \quad \xi &=0 \end{align}

and

(3.6h–j) \begin{align} f'\rightarrow 1, \quad \theta \rightarrow 1, \quad \varphi \rightarrow \phi _{\infty } \quad \textrm {as} \quad \xi \rightarrow \infty , \end{align}

where a prime denotes differentiation with respect to $\xi$ .

3.2. Boundary-layer simplifications

In the limits $Le\rightarrow \infty$ and ${\textit{Sc}}\rightarrow \infty$ , (3.6c ) simplifies to $\varphi '=0$ , implying $\varphi =\phi _{\infty }$ everywhere. Consequently, $\mu$ , $\rho$ , $c$ and $k$ are constants, and the boundary-layer equations (3.6a ) and (3.6b ) reduces to

(3.7a,b) \begin{align} 2\rho \mu f'''+ ff''=0 \quad \textrm {and} \quad 2\rho \mu \theta '' + \widehat {\textit{Pr}}f\theta ' = 0, \end{align}

where $\widehat {\textit{Pr}}=\mu cPr/k$ .

A further simplification of the boundary-layer equations is obtained by introducing

(3.8a–c) \begin{align} p=\rho \hat {p}, \quad T=1+(T_w-1)\widehat {T}, \quad \widehat {\textit{Re}} = \frac {\rho }{\mu }\textit{Re} =\frac {U_\infty ^*L^*\rho ^*}{\mu ^*}, \end{align}

into the governing (2.12) and following the procedure outlined in § 3.1 with $\widehat {\textit{Re}}\rightarrow \infty$ , to give

(3.9a,b) \begin{align} 2 f'''+ff''=0 \quad \textrm {and} \quad 2\widehat {\theta }'' + \widehat {\textit{Pr}}f\widehat {\theta }' = 0, \end{align}

subject to the boundary conditions

(3.9c–e) \begin{align} f=f'=0, \; \widehat {\theta }=1 \quad \textrm {on} \quad \xi =0, \end{align}

and

(3.9f,g) \begin{align} f'\rightarrow 1, \quad \widehat {\theta }\rightarrow 0, \quad \textrm {as} \quad \xi \rightarrow \infty . \end{align}

Here, the similarity solution $\widehat {\theta }(\xi )=\widehat {T}(x,y)$ , with $\rho ^*$ and $\mu ^*$ representing the nanofluid density (2.2) and viscosity (2.5), respectively. The rescaling in (3.8) absorbs the density and viscosity, removing them from the governing equations. Consequently, the equations simplify to the standard Blasius formulation with a modified Prandtl number $\widehat {\textit{Pr}}$ and Reynolds number $\widehat {\textit{Re}}$ based on the nanofluid quantities, allowing the nanofluid flow to be modelled as a single-phase fluid. Thus, in this simplified formulation, the flow characteristics are identical to those obtained for the standard Blasius flow, irrespective of the nanofluid quantities. Hence, in the absence of BM and TP, the Reynolds number of the nanofluid flow ${\textit{Re}}$ is given in terms of $\widehat {\textit{Re}}$ as ${\textit{Re}}=\mu \widehat {\textit{Re}}/\rho$ . A detailed description of the Navier–Stokes equations in the absence of BM and TP, leading to the derivation of (3.9), is given in Appendix A.

3.3. Boundary-layer solutions

On the left-hand side of figure 3, the steady streamwise velocity $U_{\!B}=f'(\xi )$ , temperature $T_{\!B}=\theta (\xi )$ and nanoparticle volume concentration $\phi _B=\varphi (\xi )$ are plotted for five values of $\phi _{\infty }$ and the wall temperature $T_w=2$ . Similar profiles are obtained for other values of $T_w$ . The solid, dashed and chain lines represent solutions of the full boundary-layer equations (3.6) for Cu nanoparticles in water (see table 1 for thermophysical properties). A thin concentration layer develops in the $\phi _{B}$ profile, consistent with the observations of Avramenko et al. (Reference Avramenko, Blinov and Shevchuk2011), which alters the near-wall behaviour of the velocity and temperature profiles. This behaviour is most clearly illustrated on the right-hand side of figure 3, which plots the profiles $U_{\!B}' = f''(\xi )$ , $T_{\!B}' = \theta '(\xi )$ and $\phi _B'=\varphi '(\xi )$ . These profiles reveal that, in contrast to the standard Blasius flow, $U_{\!B}'$ does not approach a constant as $\xi \rightarrow 0$ .

Figure 3.

Steady base flow profiles for variable $\phi _{\infty }$ and $T_w = 2$ , for Cu nanoparticles in water. (a) Streamwise velocity $U_{\!B}=f'(\xi )$ , (b) $U_{\!B}'=f''(\xi )$ , (c) temperature $T_{\!B}=\theta (\xi )$ , (d) $T_{\!B}'=\theta '(\xi )$ , (e) nanoparticle volume concentration $\phi _B=\varphi (\xi )$ and ( f) $\phi '_B=\varphi '(\xi )$ . Dotted lines depict the equivalent solutions in the instance $Le\rightarrow \infty$ and ${\textit{Sc}}\rightarrow \infty$ .

When BM and TP are neglected (i.e. ${\textit{Sc}}\rightarrow \infty$ and $Le\rightarrow \infty$ ), the concentration layer disappears with $\phi _B=\phi _{\infty }$ everywhere (see the vertical dotted lines in figure 3 e). In this limit, the standard Blasius flow structure is recovered, with $U_{\!B}'$ approaching a constant near the wall, as indicated by the dotted lines in figure 3(b).

Table 2 compares the base flow properties on the plate surface for varying $\phi _{\infty }$ and $T_w=2$ . The differences between the results obtained with and without BM and TP are negligible for $\phi _{\infty }\lt 10^{-3}$ , but grow, due to the impact of the concentration layer, at larger $\phi _{\infty }$ .

Table 2.

Base flow properties on $\xi =0$ for variable $\phi _{\infty }$ and $T_w=2$ , where a prime denotes differentiation with respect to the similarity variable $\xi$ . Solutions based on Cu nanoparticles in water, while the results in brackets correspond to the solutions obtained in the absence of BM and TP.

Since the base flow profiles in figure 3 are plotted against the density-weighted similarity variable $\xi$ , a physically meaningful measure of the boundary-layer thickness is provided by the displacement thickness. The dimensional displacement thickness $\delta ^*_1=x^*\delta _1/{\textit{Re}}_x^{1/2}$ and momentum thickness $\delta ^*_2=x^*\delta _2/{\textit{Re}}_x^{1/2}$ , for

(3.10a,b) \begin{align} \delta _1 = \int _0^{\infty }\frac {1}{\rho (\xi )}-\frac {f'(\xi )}{\rho _{\infty }} \;\textrm {d}\xi \quad \textrm {and} \quad \delta _2 = \int _0^{\infty }\frac {f'(\xi )}{\rho _{\infty }}\left (1-f'(\xi )\right ) \;\textrm {d}\xi , \end{align}

are shown in figure 4, along with the shape factor $H=\delta ^*_1/\delta ^*_2$ . Here, ${\textit{Re}}_x=U_{\infty }^*x^*\rho _{\textit{bf}}^*/\mu _{\textit{bf}}^*$ and $\rho _{\infty }=\rho ^*_{\infty }/\rho ^*_{\textit{bf}}$ denotes the dimensionless free stream density. Results are plotted for all seven nanoparticle materials listed in table 1. For all but two of these materials, both $\delta _1$ and $\delta _2$ decrease as $\phi _{\infty }$ increases. The most significant reductions occur for Ag and Cu nanoparticles, which have the highest densities (and the largest non-dimensional $\hat {\rho }$ values). In contrast, silicon (Si) and aluminium (Al) nanoparticles, which have the lowest densities (and the smallest values of $\hat {\rho }$ ), show an increase in $\delta _1$ and $\delta _2$ as $\phi _{\infty }$ increases. (Solutions corresponding to the case without BM and TP are nearly identical to those shown in figure 4.)

Figure 4.

(a) Displacement thickness $\delta _1$ , (b) momentum thickness $\delta _2$ and (c) shape factor $H$ as functions of the free stream nanoparticle volume concentration $\phi _{\infty }$ , for different nanoparticle materials.

The thermal displacement thickness $\delta ^*_T=x^*\delta _T/{\textit{Re}}_x^{1/2}$ and concentration displacement thickness $\delta ^*_{\phi }=x^*\delta _{\phi }/{\textit{Re}}_x^{1/2}$ , for

(3.11a,b) \begin{align} \delta _T = \int _0^{\infty }\frac {1}{\rho (\xi )}-\frac {\theta (\xi )-T_w}{\rho _{\infty }(1-T_w)} \;\textrm {d}\xi \quad \textrm {and} \quad \delta _{\phi } = \int _0^{\infty }\frac {1}{\rho (\xi )}-\frac {\varphi (\xi )}{\rho _{\infty }\phi _{\infty }} \;\textrm {d}\xi , \end{align}

are plotted in figure 5 as a function of $\phi _{\infty }$ . In contrast to the displacement thickness $\delta _1$ , the thermal displacement thickness $\delta _T$ increases with increasing $\phi _{\infty }$ for all seven nanoparticle materials. The most pronounced increases are observed for the less dense materials, Al and Si. On the other hand, the concentration displacement thickness $\delta _{\phi }$ (plotted on a semilogarithmic scale along the horizontal axis) exhibits only minor variations across the range of $\phi _{\infty }$ shown. However, noticeable differences arise between the materials. Notably, $\text{TiO}_2$ and $\text{Al}_2\text{O}_3$ exhibit larger values of $\delta _{\phi }$ than the other materials. This can be attributed to their respective $N_{{\textit{BT}}}$ values being an order of magnitude smaller than those of the other materials (see table 1). Thus, TP effects are more dominant than BM effects for these particular materials. Moreover, as $\phi _{\infty }$ approaches zero, $\delta _{\phi }$ tends towards a positive constant, indicating that $\varphi$ approaches a limiting solution. This behaviour will be examined in further detail in § 3.4.

Figure 5.

(a) Thermal displacement thickness $\delta _T$ and (b) concentration displacement thickness $\delta _{\phi }$ as functions of the free stream nanoparticle volume concentration $\phi _{\infty }$ , for different nanoparticle materials.

Figure 6.

Scaled local Nusselt number $Nu{\textit{Re}}_x^{-1/2}$ as a function of the free stream nanoparticle volume concentration $\phi _{\infty }$ , for different nanoparticle materials.

Despite the thickening of the thermal boundary-layer, the local Nusselt number, defined as

(3.12) \begin{align} {\textit{Nu}} = \frac {{\textit{Re}}_x^{1/2}\rho _w k_w\theta '(0)}{1-T_w}, \end{align}

increases with increasing $\phi _{\infty }$ , as shown in figure 6. Thus, all of the nanoparticles improve the heat transfer capabilities of the fluid. The most pronounced increases in $Nu$ are observed for denser materials with higher thermal conductivities and smaller specific heat capacities, such as Ag and Cu nanoparticles. Consequently, these materials have greater thermodynamic benefits.

3.4. Asymptotic behaviour in the limit $\phi _{\infty }\rightarrow 0$

The behaviour of the steady base flow is now examined in the limit as the free stream nanoparticle volume concentration $\phi _{\infty }$ approaches zero. Similarity variables $f$ , $\theta$ and $\varphi$ are expanded in powers of $\phi _{\infty }$ , as

(3.13a) \begin{align} f(\xi ) &= f_0(\xi ) + \phi _{\infty }f_1(\xi ) + O(\phi _{\infty }^2), \end{align}

(3.13b) \begin{align} \theta (\xi ) &= \theta _0(\xi ) + \phi _{\infty }\theta _1(\xi ) + O(\phi _{\infty }^2), \end{align}

(3.13c) \begin{align} \varphi (\xi ) &= \phi _{\infty }\varphi _1(\xi ) + O(\phi _{\infty }^2), \end{align}

while the physical quantities $\mu$ , $\rho$ , $c$ and $k$ are of the form

(3.13d) \begin{align} (\mu ,\rho ,c,k)(\xi ) = 1 +\phi _{\infty }(\mu _1,\rho _1,c_1,k_1)(\xi ) + O(\phi _{\infty }^2). \end{align}

Substituting (3.13) into (3.6a ) and (3.6b ) and retaining the leading-order terms yields the Blasius boundary-layer equations for the velocity and temperature

(3.14a,b) \begin{align} 2f_0^{\prime\prime\prime}+f_0f_0^{\prime\prime} = 0 \quad \textrm {and} \quad 2\theta _0^{\prime\prime} + \textit{Pr} f_0\theta _0^{\prime} = 0, \end{align}

subject to the boundary conditions

(3.14c–e) \begin{align}& f_0 = f^{\prime}_0 = 0, \quad \theta _0 = T_w \quad \textrm {on} \quad \xi =0, \end{align}

(3.14f,g) \begin{align}&\quad f_0^{\prime}\rightarrow 1, \quad \theta _0\rightarrow 1 \quad \textrm {as} \quad \xi \rightarrow \infty . \end{align}

Moreover, substituting (3.13) into (3.6c ) and equating terms of order $\phi _{\infty }$ gives the following second-order differential equation for $\varphi _1$ :

(3.15a) \begin{align} \theta _0\varphi _1^{\prime\prime} + \left (\theta _0^{\prime}+\frac {\theta _0^{\prime}}{N_{{\textit{BT}}}\theta _0} + \frac {Sc f_0}{2}\right )\varphi _1^{\prime} + \frac {1}{N_{{\textit{BT}}}}\left (\frac {\theta _0^{\prime\prime}}{\theta _0} - \left (\frac {\theta _0^{\prime}}{\theta _0}\right )^2 \right )\varphi _1 = 0, \end{align}

subject to the boundary conditions

(3.15b) \begin{align}& \theta _0\varphi _1^{\prime}+\frac {\varphi _1 \theta _0^{\prime}}{N_{{\textit{BT}}}\theta _0}=0 \quad \textrm {on} \quad \xi =0, \end{align}

(3.15c) \begin{align}&\qquad \varphi _1\rightarrow 1 \quad \textrm {as} \quad \xi \rightarrow \infty . \end{align}

Substituting the solution of (3.14) into (3.15) establishes the limiting boundary-value problem for $\varphi _1$ , with solutions presented in figure 7(a) for all seven nanoparticle materials given in table 1. These solutions illustrate the influence of the BM to TP ratio $N_{{\textit{BT}}}$ on the behaviour of the concentration layer. As $N_{{\textit{BT}}}$ decreases, the concentration layer becomes thicker. Notably, the solution corresponding to $\text{TiO}_2$ , represented by the green solid line, exhibits an overshoot near the wall, where $\varphi _1\gt 1$ before approaching the free stream value for larger $\xi$ (beyond the range shown in figure 7 a). Conversely, as $N_{{\textit{BT}}}$ increases and BM dominates diffusion effects, the nanoparticle volume concentration $\varphi _1\rightarrow 1$ for all $\xi$ , indicating a uniform concentration profile across the boundary layer.

Figure 7.

(a) Scaled profile of the nanoparticle volume concentration $\varphi _1$ in the limit $\phi _{\infty }\rightarrow 0$ , for different nanoparticle materials. (b,c) Comparisons between the limiting solution $\varphi _1$ and numerical solutions $\phi _B/\phi _{\infty }$ for $\phi _{\infty }=10^{-4}$ , $\phi _{\infty }=10^{-2}$ ,and $\phi _{\infty }=10^{-1}$ , for Cu and $\text{TiO}_2$ nanoparticles.

Figures 7(b) and 7(c) compare the limiting solution $\varphi _1$ and numerical solutions $\phi _B/\phi _{\infty }$ for $\phi _{\infty }\in [10^{-4},10^{-1}]$ , for Cu and $\text{TiO}_2$ nanoparticles, respectively. In both cases, the numerical solution converges to the limiting profile $\varphi _1$ as $\phi _{\infty }\rightarrow 0$ . Indeed, significant deviations only emerge for $\phi _{\infty }=10^{-1}$ .

3.5. The concentration layer

The base flow profiles in figures 3 and 7 reveal a thin concentration layer within the boundary layer, similar to the particle concentration layer reported by Pelekasis & Acrivos (Reference Pelekasis and Acrivos1995) for the flow of a well-mixed particle suspension past a flat plate. As ${\textit{Sc}}\rightarrow \infty$ , the concentration layer narrows. Since $U_{\!B}\sim Y$ as $Y\rightarrow 0$ , the following transformations are introduced to balance the diffusion and convection terms in (3.3d ):

(3.16a–c) \begin{align} Y = {\textit{Sc}}^{-1/3}\bar {Y}, \quad U_{\!B} = {\textit{Sc}}^{-1/3}\bar {U}_{\!B}, \quad V_{\!B} = {\textit{Sc}}^{-2/3}\bar {V}_{\!B}, \end{align}

which gives the rescaled concentration equation

(3.17) \begin{align} \frac {\partial (\phi _B \bar {U}_{\!B})}{\partial x} + \frac {\partial (\phi _B \bar {V}_{\!B})}{\partial \bar {Y}} = \frac {\partial }{\partial \bar {Y}}\left (T_{\!B}\frac {\partial \phi _B}{\partial \bar {Y}}+\frac {\phi _B}{N_{{\textit{BT}}}T_{\!B}}\frac {\partial T_{\!B}}{\partial \bar {Y}}\right ). \end{align}

Thus, the concentration layer has a characteristic thickness of $O({\textit{Re}}^{-1/2}{\textit{Sc}}^{-1/3})$ .

Substituting (3.16) into (3.3a )–(3.3c ), with $Le\rightarrow \infty$ and

(3.18) \begin{align} \phi _B=\phi _{\infty }+\frac {\psi (x,\bar {Y})}{{\textit{Sc}}^{1/3}}, \end{align}

gives to leading order,

(3.19a–c) \begin{align} \frac {\partial \bar {U}_{\!B}}{\partial x} + \frac {\partial \bar {V}_{\!B}}{\partial \bar {Y}} = 0, \quad \frac {\partial ^2\bar {U}_{\!B}}{\partial \bar {Y}^2} = 0, \quad \frac {\partial ^2T_{\!B}}{\partial \bar {Y}^2} = 0. \end{align}

The leading-order term in the concentration equation (3.17) is also given by (3.19a ). Thus,

(3.20a–c) \begin{align} \bar {U}_{\!B} = \frac {\hat {\lambda }\bar {Y}}{x^{1/2}}, \quad \bar {V}_{\!B} = \frac {\hat {\lambda }\bar {Y}^2}{4x^{3/2}}, \quad T_{\!B} = T_w + \frac {\hat {\sigma }\bar {Y}}{{\textit{Sc}}^{1/3}x^{1/2}}, \end{align}

for $\hat {\lambda }=\rho _{w}f''(0)$ and $\hat {\sigma }=\rho _{w}\theta '(0)$ .

The next order term in the concentration equation (3.17) is given as

(3.21a) \begin{align} \bar {U}_{\!B}\frac {\partial \psi }{\partial x} + \bar {V}_{\!B}\frac {\partial \psi }{\partial \bar {Y}} = T_w\frac {\partial ^2\psi }{\partial \bar {Y}^2}, \end{align}

with boundary conditions

(3.21b) \begin{align} \frac {\partial \psi }{\partial \bar {Y}} +\frac {\phi _{\infty }\hat {\sigma }}{N_{{\textit{BT}}}T_w^2x^{1/2}}=0 \quad \textrm {on} \quad \bar {Y}=0 \end{align}

and

(3.21c) \begin{align} \psi \rightarrow 0 \quad \textrm {as} \quad \bar {Y}\rightarrow \infty . \end{align}

Introducing the similarity transformation

(3.22a) \begin{align} \psi (x,\bar {Y}) =\frac {\phi _{\infty }\hat {\sigma }\varPsi (\bar {\eta })}{N_{{\textit{BT}}}\hat {\lambda }^{1/3}T_w^{5/3}}, \end{align}

for

(3.22b) \begin{align} \bar {\eta } = \left (\frac {\hat {\lambda }}{T_w}\right )^{1/3}\frac {\bar {Y}}{x^{1/2}}, \end{align}

gives the similarity equation

(3.23a) \begin{align} \frac {{d}^2\varPsi }{\textrm {d}\bar {\eta }^2}+\frac {\bar {\eta }^2}{4}\frac {\textrm {d}\varPsi }{\textrm {d}\bar {\eta }}=0, \end{align}

with boundary conditions

(3.23b) \begin{align} \frac {\textrm {d}\varPsi }{\textrm {d}\bar {\eta }}=-1 \quad \textrm {on} \quad \bar {\eta }=0 \end{align}

and

(3.23c) \begin{align} \varPsi \rightarrow 0 \quad \textrm {as} \quad \bar {\eta }\rightarrow \infty . \end{align}

The solution for $\varPsi$ is given in terms of the upper incomplete Gamma function $\varGamma$ ,

(3.24) \begin{align} \varPsi (\bar {\eta }) = \left (\frac {2}{3}\right )^{2/3}\varGamma \left (\frac {1}{3},\frac {\bar {\eta }^3}{12}\right )\!, \end{align}

and is plotted in figure 8(a). At the wall, $\varPsi (0) \approx 2.0444$ . Hence, to a first approximation, the nanoparticle volume concentration is given by

(3.25) \begin{align} \phi _B = \phi _{\infty }\left (1+\frac {\hat {\sigma }\varPsi (\bar {\eta })}{N_{{\textit{BT}}}\hat {\lambda }^{1/3}T_w^{5/3}{\textit{Sc}}^{1/3}} \right ). \end{align}

Figure 8.

(a) Similarity solution $\varPsi$ for the nanoparticle volume concentration, as given by (3.24). (b–e) Nanoparticle volume concentration profiles $\phi _B$ given by the exact solution to (3.6) (solid blue lines) and the approximate solution (3.25) (dashed red), for Cu nanoparticles.

Figures 8(b) and 8(c) compare the exact nanoparticle volume concentration profiles $\phi _B$ , obtained by solving (3.6), with the approximate solution given by (3.25), for Cu nanoparticles and $T_w=2$ . Results are plotted for $\phi _{\infty }=10^{-3}$ and $\phi _{\infty }=10^{-2}$ . In both cases, the approximate solution is qualitatively similar to the exact solution, with only minor differences near the wall, corresponding to a maximum relative error of approximately 3 %. Such small differences are to be expected since $N_{{\textit{BT}}}Sc^{1/3}\sim O(1)$ for the parameter settings used in figures 8(b) and 8(c). For materials with smaller $N_{{\textit{BT}}}$ values, such as $\text{Al}_2\text{O}_3$ and $\text{TiO}_2$ , the approximation is less accurate, and higher-order terms are required to improve the solution. However, by increasing both ${\textit{Sc}}$ and $N_{{\textit{BT}}}$ , as is modelled in figures 8(d) and 8(e), the agreement between the exact and approximate solutions improves significantly, with the maximum relative error reduced to 0.001 %.

4.1. Linearised stability equations

The linear stability equations are derived by decomposing the total velocity, pressure, temperature and nanoparticle volume concentration fields as

(4.1a–f) \begin{align} u &{}={} U_{\!B} + \epsilon \tilde {u}, \quad v{}={} {\textit{Re}}^{-1/2}V_{\!B} + \epsilon \tilde {v}, \quad w{}={} \epsilon \tilde {w}, \nonumber \\ p &{}={} \epsilon \tilde {p}, \qquad \quad\,\, T{}={} T_{\!B}+\epsilon \tilde {T}, \qquad \quad\,\, \phi {}={} \phi _B+\epsilon \tilde {\phi }, \end{align}

for perturbations $\tilde {\boldsymbol{q}}=(\tilde {\boldsymbol{u}}, \tilde {p}, \tilde {T}, \tilde {\phi })$ , with $\tilde {\boldsymbol{u}}=(\tilde {u}, \tilde {v}, \tilde {w})$ and $\epsilon \ll 1$ . Similarly,

(4.2a–e) \begin{align} \rho &= \rho _B+\epsilon \tilde {\rho }, \quad \rho c {}={} (\rho c)_B+\epsilon \tilde {\rho }\tilde {c}, \quad c {}={} c_{\!B}+\epsilon \tilde {c}, \nonumber \\ \mu &{}={} \mu _B+\epsilon \tilde {\mu }, \quad k {}={} k_B+\epsilon \tilde {k}. \end{align}

Here, base flow quantities $\boldsymbol{Q}_{\!B}=(U_{\!B}, V_{\!B}, T_{\!B}, \phi _B)$ depend on $x$ and $y$ , while perturbations $\tilde {\boldsymbol{q}}$ are functions of $\boldsymbol{x}$ and $t$ . Substituting (4.1) and (4.2) into (2.12), and linearising in $\epsilon$ , gives the following linear stability equations:

(4.3a) \begin{align}&\qquad\qquad\qquad \rho _B\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}} + \frac {\partial \tilde {\rho }}{\partial t}+U_{\!B}\frac {\partial \tilde {\rho }}{\partial x}+\rho _{B,y}\tilde {v} = g_1(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

(4.3b) \begin{align}&\qquad \rho _B\left (\frac {\partial \tilde {u}}{\partial t}+U_{\!B}\frac {\partial \tilde {u}}{\partial x}+U_{B,y}\tilde {v}\right ) = - \frac {\partial \tilde {p}}{\partial x} + \frac {1}{\textit{Re}}\left (\mu _B\bigg ({\nabla} ^2\tilde {u} + \frac {1}{3}\frac {\partial }{\partial x}\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}\right ) \nonumber\\&\qquad\quad + \mu _{B,y}\left (\frac {\partial \tilde {v}}{\partial x} + \frac {\partial \tilde {u}}{\partial y}\right ) + U_{B,yy}\tilde {\mu } + U_{B,y}\frac {\partial \tilde {\mu }}{\partial y}\Bigg ) + g_2(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

(4.3c) \begin{align}& \rho _B\left (\frac {\partial \tilde {v}}{\partial t}+U_{\!B}\frac {\partial \tilde {v}}{\partial x}\right ) = - \frac {\partial \tilde {p}}{\partial y} + \frac {1}{\textit{Re}}\left (\mu _B\bigg ({\nabla} ^2\tilde {v} + \frac {1}{3}\frac {\partial }{\partial y}\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}\right ) \nonumber \\&\qquad\qquad\qquad\qquad\quad + \frac {2\mu _{B,y}}{3}\left (2\frac {\partial \tilde {v}}{\partial y} - \left (\frac {\partial \tilde {u}}{\partial x}+\frac {\partial \tilde {w}}{\partial z}\right )\right ) + U_{B,y}\frac {\partial \tilde {\mu }}{\partial x}\Bigg ) + g_3(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

(4.3d) \begin{align}&\qquad\qquad \rho _B\left (\frac {\partial \tilde {w}}{\partial t}+U_{\!B}\frac {\partial \tilde {w}}{\partial x}\right ) = - \frac {\partial \tilde {p}}{\partial z} + \frac {1}{\textit{Re}}\Bigg (\mu _B\left ({\nabla} ^2\tilde {w} + \frac {1}{3}\frac {\partial }{\partial z}\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}\right )\nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\quad + \mu _{B,y}\left (\frac {\partial \tilde {v}}{\partial z} + \frac {\partial \tilde {w}}{\partial y}\right )\Bigg ) + g_4(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

(4.3e) \begin{align}&\quad \rho _BT_{\!B}\left (\frac {\partial \tilde {c}}{\partial t}+U_{\!B}\frac {\partial \tilde {c}}{\partial x}+c_{B,y}\tilde {v}\right ) + (\rho c)_B\left (\frac {\partial \tilde {T}}{\partial t}+U_{\!B}\frac {\partial \tilde {T}}{\partial x}+T_{B,y}\tilde {v}\right ) \nonumber \\&\qquad = \frac {1}{\textit{Re} \textit{Pr}}\left (\frac {\partial }{\partial y}\left (k_B\frac {\partial \tilde {T}}{\partial y}+T_{B,y}\tilde {k}\right ) + k_B\widehat {{\nabla} }^2\tilde {T}\right ) \nonumber \\&\qquad\quad + \frac {1}{\textit{Re} \textit{Pr} \textit{Le}}\left (T_{B,y}\mathcal{A} + \mathcal{B}\frac {\partial \tilde {T}}{\partial y}\right ) + g_5(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

(4.3f) \begin{align}&\quad \phi _B\boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}+\frac {\partial \tilde {\phi }}{\partial t}+U_{\!B}\frac {\partial \tilde {\phi }}{\partial x}+\phi _{B,y}\tilde {v} = \frac {1}{\textit{Re} \textit{Sc}}\left (\frac {\partial \mathcal{A}}{\partial y}+T_{\!B}\widehat {{\nabla} }^2\tilde {\phi } + \frac {\phi _B}{N_{{\textit{BT}}}T_{\!B}}\hat {{\nabla} }^2\tilde {T} \right )\nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad + g_6(V_{\!B},\boldsymbol{Q}_{B,x}), \end{align}

where functions $g_{\star }$ depend on the wall-normal velocity $V_{\!B}$ and $x$ -derivatives of the base flow $\boldsymbol{Q}_{\!B}$ , and

(4.4) \begin{align}& \mathcal{A} = T_{B}\frac {\partial \tilde {\phi }}{\partial y}+\phi _{B,y}\tilde {T}+\frac {1}{N_{{\textit{BT}}}T_{\!B}}\left (\phi _B\frac {\partial \tilde {T}}{\partial y}+T_{B,y}\tilde {\phi }-\frac {\phi _BT_{B,y}}{T_{\!B}}\tilde {T}\right )\!, \end{align}

(4.5) \begin{align}& \mathcal{B} = \phi _{B,y}T_{\!B}+\frac {\phi _BT_{B,y}}{N_{{\textit{BT}}}T_{\!B}} \end{align}

and

(4.6) \begin{align} \widehat {{\nabla} }^2=\frac {\partial ^2}{\partial x^2}+\frac {\partial ^2}{\partial z^2}. \end{align}

(The exact form of the functions $g_{\star }$ are given in Appendix B.) The corresponding boundary conditions are given as

(4.7a–e) \begin{align} \tilde {u}=\tilde {v}=\tilde {w}=\tilde {T}=\mathcal{A}=0 \quad \textrm {on} \quad y=0, \end{align}

and

(4.7f–k) \begin{align} \tilde {u}\rightarrow 0, \; \tilde {v}\rightarrow 0, \; \tilde {w}\rightarrow 0, \; \tilde {p}\rightarrow 0, \; \tilde {T}\rightarrow 0, \; \tilde {\phi }\rightarrow 0 \quad \textrm {as} \quad y\rightarrow \infty . \end{align}

The length scale $L^*$ used in the subsequent linear stability analysis is based on the displacement thickness $\delta _1^*$ , to give the Reynolds number

(4.8) \begin{align} {\textit{R}}=\frac {U_{\infty }^*\delta _1^*\rho _{\textit{bf}}^*}{\mu _{\textit{bf}}^*}, \end{align}

which ensures consistency with earlier investigations (Mack Reference Mack1984; Schmid & Henningson Reference Schmid and Henningson2001). This gives the following relationships: ${\textit{R}} = \delta _1{\textit{Re}}_x^{1/2}$ and ${\textit{R}} = \delta _1(xRe)^{1/2}$ . Consequently, $Re$ in the system of equations (4.3) is replaced with $R$ .

Additionally, the parallel flow approximation is imposed, where the flow is assumed to be in the $x$ -direction and depends only on the wall-normal $y$ -direction, i.e. $g_{\star }=0$ . Subsequently, perturbations $\tilde {\boldsymbol{q}}$ are decomposed into the normal mode form

(4.9) \begin{align} \tilde {\boldsymbol{q}}(\boldsymbol{x},t) = \breve {\boldsymbol{q}}(y)\exp {(\textrm {i}(\alpha x+\beta z-\omega t))} + \textrm {c.c.}, \end{align}

(and similarly for quantities $\tilde {\rho }$ , $\tilde {\mu }$ , etc.) for a streamwise wavenumber $\alpha \in \mathbb{R}$ , spanwise wavenumber $\beta \in \mathbb{R}$ , and frequency $\omega \in \mathbb{C}$ . Here, $\textrm {c.c.}$ denotes the complex conjugate. Consequently, (4.3) become

(4.10a) \begin{align}&\qquad\qquad\qquad \rho _B\left (\textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right )+\textrm {D}\breve {v} \right ) + \textrm {i}\left (\alpha U_{\!B}-\omega \right )\breve {\rho }+\rho _{B,y}\breve {v} = 0, \end{align}

(4.10b) \begin{align}&\qquad\quad \rho _B\left (\textrm {i}\left (\alpha U_{\!B}-\omega \right )\breve {u}+U_{B,y}\breve {v}\right ) = -\textrm {i}\alpha \breve {p} + \frac {1}{R}\Bigg (\mu _B\bigg (\left (\textrm {D}^2 - \left (\alpha ^2+\beta ^2\right )\right )\breve {u} \nonumber \\&\qquad\qquad + \frac {\textrm {i}\alpha }{3}\left (\textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right )+\textrm {D}\breve {v}\right )\bigg ) + \mu _{B,y}\left (\textrm {i}\alpha \breve {v} + \textrm {D}\breve {u}\right ) + \left (U_{B,yy}+U_{B,y}\textrm {D}\right )\breve {\mu }\Bigg ), \end{align}

(4.10c) \begin{align}& \textrm {i}\rho _B\left (\alpha U_{\!B}-\omega \right )\breve {v} = - \textrm {D}\breve {p} + \frac {1}{R}\Bigg (\mu _B\bigg (\left (\textrm {D}^2 - \left (\alpha ^2+\beta ^2\right )\right )\breve {v} + \frac {\textrm {D}}{3}\left (\textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right ) + \textrm {D}\breve {v}\right )\bigg ) \nonumber \\&\qquad\qquad\qquad\qquad + \frac {2\mu _{B,y}}{3}\left (2\textrm {D}\breve {v} - \textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right )\right ) + \textrm {i}\alpha U_{B,y}\breve {\mu }\Bigg ), \end{align}

(4.10d) \begin{align}&\qquad \textrm {i}\rho _B\left (\alpha U_{\!B}-\omega \right )\breve {w} = -\textrm {i}\beta \breve {p} + \frac {1}{R}\Bigg (\mu _B\bigg (\left (\textrm {D}^2-\left (\alpha ^2+\beta ^2\right )\right )\breve {w} \nonumber \\&\qquad\qquad\qquad\qquad\qquad + \frac {\textrm {i}\beta }{3}\left (\textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right )+\textrm {D}\breve {v}\right )\bigg ) + \mu _{B,y}\left (\textrm {i}\beta \breve {v} + \textrm {D}\breve {w}\right )\Bigg ), \end{align}

(4.10e) \begin{align}&\quad \rho _BT_{\!B}\left (\textrm {i}\left (\alpha U_{\!B}-\omega \right )\breve {c}+c_{B,y}\breve {v}\right ) + (\rho c)_B\left (\textrm {i}\left (\alpha U_{\!B}-\omega \right )\breve {T}+T_{B,y}\breve {v}\right ) \nonumber \\&\qquad = \frac {1}{R \textit{Pr}}\left (\textrm {D}\left (k_B\textrm {D}\breve {T}+T_{B,y}\breve {k}\right ) - \left (\alpha ^2+\beta ^2\right ) k_B\breve {T}\right ) + \frac {1}{R \textit{Pr} \textit{Le}}\left (T_{B,y}\mathcal{A} + \mathcal{B}\textrm {D}\breve {T}\right )\!, \end{align}

(4.10f) \begin{align}&\qquad\qquad\qquad \phi _B\left (\textrm {i}\left (\alpha \breve {u}+\beta \breve {w}\right )+\textrm {D}\breve {v}\right )+\textrm {i}\left (\alpha U_{\!B}-\omega \right )\breve {\phi }+\phi _{B,y}\breve {v}\nonumber \\&\qquad\qquad\qquad\quad = \frac {1}{R \textit{Sc}}\left (\textrm {D}\mathcal{A}-\left (\alpha ^2+\beta ^2\right )\left (T_{\!B}\breve {\phi } + \frac {\phi _B}{N_{{\textit{BT}}}T_{\!B}}\breve {T}\right ) \right )\!, \end{align}

where $\textrm {D}={d}/\textrm {d}y$ . The exact form of the perturbed quantities, including $\breve {\rho }$ , $\breve {\mu }$ , etc., are given in Appendix B.

4.2. Numerical methods

A temporal linear stability analysis was conducted using the Chebyshev collocation method developed by Trefethen (Reference Trefethen2000). Derivatives in the $y$ -direction were approximated using Chebyshev matrices, with $N$ Chebyshev mesh points mapped from the semi-infinite physical domain $y \in [0, \infty )$ onto the computational interval $\zeta \in [1,-1]$ via the coordinate transformation

(4.11) \begin{align} y = \frac {l(1-\zeta )}{1+\zeta }, \end{align}

where $l$ is a stretching parameter.

The linear stability equations (4.10) were transformed into the following eigenvalue problem:

(4.12) \begin{align} \unicode{x1D63C}\breve {\boldsymbol{q}}^T= \omega \unicode{x1D63D}\breve {\boldsymbol{q}}^T, \end{align}

where $\unicode{x1D63C}$ and $\unicode{x1D63D}$ are $6N\times 6N$ matrices. The frequencies $\omega$ and the corresponding linear perturbations $\breve {\boldsymbol{q}}$ were then computed using the eig command in MATLAB.

Table 3 presents the frequency $\omega$ corresponding to the TS wave for varying values of $N$ and $l$ , for Cu nanoparticles and free stream nanoparticle volume concentrations $\phi _{\infty }\in [10^{-4},10^{-1}]$ . In each case, the Reynolds number ${\textit{R}}=500$ , the streamwise wavenumber $\alpha =0.3$ , the spanwise wavenumber $\beta =0$ and the wall temperature $T_w=2$ . The results are identical to four decimal places for all $l$ considered when $N\geqslant 64$ , indicating that the thin concentration layer is well-resolved and the frequencies $\omega$ have converged. Therefore, for the remainder of this investigation, $N=96$ Chebyshev mesh points were used with the mapping parameter $l=2$ .

Table 3.

Frequencies $\omega =\omega _r+\textrm {i}\omega _i$ for variable $N$ and $l$ , for ${\textit{R}}=500$ , $\alpha =0.3$ , $\beta =0$ , $T_w=2$ and $\phi _{\infty }=10^{-4}$ , $\phi _{\infty }=10^{-2}$ , $\phi _{\infty }=10^{-1}$ . Here, $\omega _i\gt 0$ corresponds to linearly unstable behaviour.

4.3. Numerical results

In the following linear stability analysis, unless stated otherwise, the nanofluid is composed of Cu nanoparticles dispersed in a base fluid of water. In addition, the wall temperature $T_w=2$ . (The case $T_w=2$ was selected as a representative case. However, as shown in Appendix C, variations in wall temperature have negligible influence on the results.)

4.3.1. Eigenspectrum

Figure 9 presents a representative eigenspectrum in the complex $\omega$ -plane for the parameter settings $R=500$ , $\alpha =0.3$ and $\beta =0$ , and three values of $\phi _{\infty }$ . For the standard Blasius flow without nanoparticles, these conditions are linearly stable. The left-hand plots display the eigenspectrum on a large scale, while the right-hand plots provide a zoomed-in view. The blue circular markers correspond to solutions where BM and TP are ignored, whereas the red crosses indicate the corresponding solutions when these effects are included. The black star markers represent the eigenspectrum for the Blasius flow without nanoparticles, where the nanoparticle volume concentration equations have been removed from the analysis.

Figure 9.

Eigenspectrum in the $(\omega _r,\omega _i)$ -plane for $R=500$ , $\alpha =0.3$ , $\beta =0$ , $T_w=2$ , and (a,b) $\phi _{\infty }=10^{-4}$ , (c,d) $\phi _{\infty }=10^{-3}$ and (e, f) $\phi _{\infty }=10^{-2}$ . Black asterisk markers represent solutions of the Blasius flow, while blue circles and red crosses represent solutions of the nanofluid flow without (BM/TP off) and with (BM/TP on) BM and TP.

Consistent with previous studies (Mack Reference Mack1976; Grosch & Salwen Reference Grosch and Salwen1978; Salwen & Grosch Reference Salwen and Grosch1981; Schmid & Henningson Reference Schmid and Henningson2001), the eigenspectrum consists of multiple branches. A discrete set of modes are located on the A-branch (Mack Reference Mack1976) in the upper left-hand corner of figure 9(a,c,e). This branch contains the TS wave, which is highlighted in the right-hand plots and discussed further below. Additionally, the eigenspectrum features three continuous branches, each associated with different governing equations. (The eigenspectrum shown is a discrete representation of the continuous spectrum, with the resolution governed by the number of Chebyshev mesh points $N$ .) The first two branches, approximately aligned with the vertical axis, are associated with the momentum and energy equations, respectively. As the number of Chebyshev mesh points $N$ increases, these two branches shift to the right towards the vertical line $\omega _r\rightarrow \alpha$ , although their qualitative behaviour is unchanged. The third continuous branch, associated with the nanoparticle volume concentration equation, runs parallel to the real $\omega$ -line but with a negative imaginary part. Like the other two continuous branches, this branch also shifts to the right as $N$ increases, but at a significantly slower rate due to the size of the Schmidt number ${\textit{Sc}}$ . Notably, when BM and TP are neglected, this branch is located along the real $\omega$ -line (i.e. $\omega _i=0$ ), as expected, since equation (4.10f ) simplifies to

(4.13) \begin{align} \left (\alpha U_{\!B}-\omega \right )\breve {\phi } = 0 \end{align}

in this case.

The zoomed-in plots on the right-hand side of figure 9 focus on the behaviour of the frequency $\omega$ of the TS wave as the free stream nanoparticle volume concentration $\phi _{\infty }$ increases. For $\phi _{\infty }=10^{-4}$ , the value of $\omega$ closely matches that of the Blasius flow without nanoparticles, with linearly stable conditions, as the imaginary part of $\omega$ is negative. However, as $\phi _{\infty }$ increases, a noticeable shift occurs. At $\phi _{\infty }=10^{-3}$ , the frequency $\omega$ shifts slightly to the left and upward in the $\omega$ -plane, remaining linearly stable but less stable than the standard Blasius flow. With a further increase to $\phi _{\infty }=10^{-2}$ , $\omega$ moves into the upper half-plane, where a positive imaginary part indicates linearly unstable behaviour. Thus, for the given flow conditions, the nanofluid destabilises the TS wave. Furthermore, the differences in $\omega$ obtained with and without the effects of BM and TP are minimal, with only slight variations in the real component and no discernible changes in the imaginary component. (In addition to the frequency $\omega$ of the TS wave, eigenspectra from the branch arising from the nanoparticle volume concentration equation are also shown in figure 9(b,d), further illustrating how this branch aligns with the real $\omega$ -axis.)

Figure 10.

Frequency $\omega =\omega _r+\textrm {i}\omega _i$ as a function of $\phi _{\infty }$ for $R=500$ , $\alpha =0.3$ , $\beta =0$ and $T_w=2$ . (a) Real part and (b) imaginary part. The solid blue and dashed red lines represent solutions of the nanofluid flow without (BM/TP off) and with (BM/TP on) BM and TP. The horizontal chain lines indicate the corresponding solutions for the Blasius flow without nanoparticles.

Figure 10 further illustrates the variation of the frequency $\omega$ of the TS wave as the free stream nanoparticle volume concentration $\phi _{\infty }$ increases, for the same conditions as given in figure 9. The plots show the evolution of both the real and imaginary components of $\omega$ with increasing $\phi _{\infty }$ , supporting the trend observed in figure 9. As more nanoparticles are added to the base fluid, the TS wave becomes increasingly destabilised, with the imaginary part of $\omega$ shifting from negative to positive values near $\phi _{\infty }=0.008$ , signalling the onset of linear instability. Additionally, solutions demonstrate that the effects of BM and TP are negligible, since the differences between cases without (solid blue lines) and with (dashed red) these effects are minimal, with only slight variations in the real part of $\omega$ and no significant impact on the imaginary part.

4.3.2. Three-dimensional instabilities

Figure 11.

Temporal growth rate $\omega _i$ as a function of the streamwise wavenumber $\alpha$ for $R=600$ , $T_w=2$ , $\beta \in [0,0.1]$ and (a) $\phi _{\infty }=10^{-4}$ , (b) $\phi _{\infty }=10^{-3}$ and (c) $\phi _{\infty }=10^{-2}$ .

Although Squire’s theorem cannot be applied directly to the full linear stability equations (4.10), it is applicable to the simplified linear stability equations that neglect BM and TP. Since these diffusion effects have a minimal impact on both the base flow and the linear stability calculations, we conclude that Squire’s theorem is approximately valid for the full equations. Consequently, it is sufficient to limit the stability analysis to two-dimensional instabilities.

This conclusion is supported by the results shown in figure 11, which plots the temporal growth rate $\omega _i$ as a function of the streamwise wavenumber $\alpha$ , for the Reynolds number ${\textit{R}}=600$ , spanwise wavenumbers $\beta \in [0,0.1]$ and nanoparticle volume concentrations $\phi _{\infty }\in [10^{-4},10^{-2}]$ . The results indicate that $\omega _i$ decreases as $\beta$ increases, confirming that two-dimensional instabilities are more unstable than three-dimensional instabilities. Therefore, based on this and further observations, the remainder of this study focuses on two-dimensional disturbances by setting $\beta = 0$ .

4.3.3. Conditions for neutral stability

Figure 12.

Neutral stability curves in the $(R,\omega )$ -plane for variable $\phi _{\infty }$ , $\beta =0$ , $T_w=2$ and (a) Cu nanoparticles and (b) Al nanoparticles.

The neutral conditions $(\omega , R)$ for linear instability were computed using streamwise wavenumber increments of $\Delta \alpha =10^{-4}$ . To accurately trace the frequency $\omega$ associated with the TS wave within the complex $\omega$ -plane, small Reynolds number steps $\Delta {\textit{R}}=0.01$ were used. This ensured that the TS frequency was correctly identified, minimising interference with the eigenspectra found on the branch due to the nanoparticle volume concentration equation. The critical Reynolds number for the Blasius flow, in the absence of nanoparticles, was obtained as ${\textit{R}}_c\approx 519.4$ for a streamwise wavenumber $\alpha _c\approx 0.304$ , frequency $\omega _c\approx 0.121$ and phase speed $s_c=\omega _c/\alpha _c\approx 0.397$ , in agreement with previous studies (Schmid & Henningson Reference Schmid and Henningson2001).

Neutral stability curves were obtained for free stream nanoparticle volume concentrations $\phi _{\infty }\in [0,4\times 10^{-2}]$ , with solutions for the Cu nanoparticles shown in figure 12(a). The destabilisation of the TS wave is further demonstrated, with neutral stability curves shifting horizontally to the left and smaller Reynolds numbers as $\phi _{\infty }$ increases. Notably, there is no discernible vertical variation in the neutral stability curves. Thus, while the critical Reynolds number $\textit{R}_c$ shrinks, the corresponding frequency $\omega _c$ , the streamwise wavenumber $\alpha _c$ and the phase velocity $s_c$ , remain relatively constant for the range of $\phi _{\infty }$ considered.

A second set of neutral stability curves is shown in figure 12(b), but for nanoparticles made of Al. Like the Cu nanoparticles, there is no vertical variation as $\phi _{\infty }$ increases. However, a small stabilising effect is observed, with neutral curves shifting to the right and marginally larger Reynolds numbers ${\textit{R}}$ . Therefore, the type of material used for the nanoparticles plays a significant role in determining whether the TS wave is stabilised or destabilised.

Figure 13 presents further evidence of the stabilising benefits of Al nanoparticles compared with the destabilising effects of Cu nanoparticles. The circular (Cu) and diamond (Al) markers indicate the critical Reynolds numbers ${\textit{R}}_c$ obtained from the full linear stability equations (4.10), with a noticeable reduction in ${\textit{R}}_c$ for Cu nanoparticles and a small increase for Al nanoparticles. Additionally, the critical Reynolds number ${\textit{R}}_c$ for these two types of nanoparticles is plotted when BM and TP are neglected, as represented by the solid blue and dashed red curves. In this case, the critical Reynolds number ${\textit{R}}_c=\mu \widehat {\textit{R}}_c/\rho$ , where $\widehat {R}_c\approx 519.4$ is the critical Reynolds number for the Blasius flow without nanoparticles. Thus, using the definition for density $\rho$ and the Brinkman dynamic viscosity $\mu$ , given by (2.13a ) and (2.14), respectively, the critical Reynolds for the nanofluid flow is approximated as

(4.14) \begin{align} {\textit{R}}_c = \frac {519.4}{(1-\phi _{\infty })^{2.5}(1+(\hat {\rho }-1)\phi _{\infty })}. \end{align}

Unsurprisingly, the results with and without BM and TP are nearly identical. Thus, the impact of these diffusion effects on the linear stability of the nanofluid flow are negligible. Table 4 lists critical Reynolds numbers ${\textit{R}}_c$ at select $\phi _{\infty }$ values for both Cu and Al nanoparticles.

Table 4.

Critical Reynolds numbers ${\textit{R}}_c$ for Cu and Al nanoparticles in a base fluid of water, while the results in brackets correspond to the solutions obtained in the absence of BM and TP.

Figure 13.

Critical Reynolds number ${\textit{R}}_c$ as a function of $\phi _{\infty }$ , for Cu nanoparticles (solid blue line and circular markers) and Al nanoparticles (dashed red line and diamond markers) in a base fluid of water without (BM/TP off) and with (BM/TP on) BM and TP.

Figure 14.

Plots of the critical Reynolds number ${\textit{R}}_c$ for the seven nanoparticle materials tabulated in table 1 in a base fluid of water, with the dynamic viscosity $\mu$ based on the Brinkman (Reference Brinkman1952) model (2.14). ( $a$ ) Here ${\textit{R}}_c$ as a function of $\phi _{\infty }$ . ( $b$ ) Contours of ${\textit{R}}_c$ in the ( $\phi _{\infty },\hat {\rho }$ )-plane, where the solid red contour represents the contour level ${\textit{R}}_c=519.4$ , matched to the critical conditions for the Blasius flow without nanoparticles.

Consequently, the critical Reynolds number $\textit{R}_c$ is governed by the dynamic viscosity $\mu$ and the density $\rho$ of the nanofluid, which are in turn influenced by the free stream nanofluid volume concentration $\phi _{\infty }$ and the ratio of densities $\hat {\rho }$ . Figure 14 illustrates $\textit{R}_c$ as approximated by equation (4.14). In the first plot, figure 14(a), $\textit{R}_c$ is plotted as a function of $\phi _{\infty }$ and demonstrates the influence of both $\phi _{\infty }$ and the material used for the nanoparticles. Denser materials with larger $\hat {\rho }$ ratios, like Ag and Cu, have a destabilising effect, while lighter materials, like Si and Al, stabilise the flow. On the other hand, $\text{Al}_2\text{O}_3$ exhibits a marginally destabilising effect at small $\phi _{\infty }$ , with a stabilising benefit realised for large $\phi _{\infty }$ (for $\phi _{\infty }\gtrapprox 0.09$ ).

Figure 14(b) further demonstrates the impact of nanofluids on the onset of linear instability, with $\textit{R}_c$ plotted in the $(\phi _{\infty },\hat {\rho })$ -plane. The solid red contour corresponds to $R_c=519.4$ (i.e. the onset of linear instability in the standard Blasius flow), with solutions illustrating the negative impact of most nanoparticle materials, except Si and Al, on the hydrodynamic stability of the flow. More specifically, for a base fluid of water, only nanoparticles with a density ratio $\hat {\rho }\lessapprox 3.5$ are stabilising.

5.1. The main deck

Here $y=\varepsilon ^4 y_2$ , for $y_2=O(1)$ , where perturbations $\tilde {\boldsymbol{q}}=(\tilde {u},\tilde {v},\tilde {p},\tilde {T},\tilde {\phi })$ are expanded as

(5.3a–e) \begin{align} \tilde {u} &{}={} \left (u_2 + O(\varepsilon )\right )E, \qquad \quad \tilde {v} {}={} \left (\varepsilon v_2+O\big(\varepsilon ^2\big)\right )E, \nonumber \\ \tilde {p} &{}={} \left (\varepsilon p_2+O\big(\varepsilon ^2\big)\right )E, \quad \tilde {T} {}={} \left (T_2+O(\varepsilon )\right )E, \nonumber \\ \tilde {\phi } &{}={} \left (\phi _2+O(\varepsilon )\right )E, \end{align}

where $u_2=u_2(x,y_2)$ , etc. Similar expansions are given for the perturbed quantities $\tilde {\mu }$ , $\tilde {\rho }$ , $\tilde {c}$ and $\tilde {k}$ . In addition, the nanoparticle volume concentration $\phi _B\sim \phi _{\infty }$ .

Substituting (5.3) into the linear stability equations (4.3) and collecting the leading-order terms, gives the solution

(5.4a–c) \begin{align} u_2 = A(x)U_{B,y_2}, \quad v_2 = -\textrm {i}\alpha _1 A(x) U_{\!B} \quad \textrm {and} \quad p_2 = p_2(x), \end{align}

where $p_2(x)$ and $A(x)$ are unknown, slowly varying, amplitude functions, representing pressure and negative displacement perturbations, respectively. Similarly,

(5.4d,e) \begin{align} T_2 = A(x)T_{B,y_2} \quad \textrm {and} \quad \phi _2 = 0. \end{align}

5.2. The upper deck

Here $y=\varepsilon ^3 y_1$ , for $y_1=O(1)$ . To match with the main deck, perturbations are expanded as

(5.5a–e) \begin{align} \tilde {u} &{}={} \left (\varepsilon u_1+O\big(\varepsilon ^2\big)\right )E, \quad \tilde {v} {}={} \left (\varepsilon v_1+O\big(\varepsilon ^2\big)\right )E, \nonumber \\ \tilde {p} &{}={} \left (\varepsilon p_1+O\big(\varepsilon ^2\big)\right )E, \quad \tilde {T} {}={} \left (\varepsilon T_1+O\big(\varepsilon ^2\big)\right )E, \nonumber \\ \tilde {\phi } &{}={} \left (\varepsilon \phi _1+O\big(\varepsilon ^2\big)\right )E, \end{align}

where $u_1=u_1(x,y_1)$ , etc. Similar expansions are again given for the perturbed quantities $\tilde {\mu }$ , $\tilde {\rho }$ , $\tilde {c}$ and $\tilde {k}$ . In addition, the base flow is effectively given by the uniform free stream conditions

(5.6a–f) \begin{align} U_{\!B} &{}\approx {} 1, \qquad V_{\!B} {}\approx {} 0, \qquad T_{\!B} {}\approx {} 1, \nonumber \\ \phi _B &{}\approx {} \phi _{\infty }, \quad\, c_{\!B} {}\approx {} c_{\infty }, \quad\, \rho _B {}\approx {} \rho _{\infty }. \end{align}

Substituting (5.5) and (5.6) into the linear stability equations (4.3), gives

(5.7) \begin{align} \left (\frac {\partial ^2}{\partial y_1}-\alpha _1^2\right )p_1=0, \end{align}

with the bounded solution as $y_1\rightarrow \infty$ given by

(5.8a) \begin{align} p_1 = P_1(x) \textrm {e}^{-\alpha _1 y_1}, \end{align}

where $P_1(x)$ is an unknown function of $x$ and $\alpha _1\gt 0$ . Moreover,

(5.8b–e) \begin{align} u_1 = -\frac {P_1(x)\textrm {e}^{-\alpha _1 y_1}}{\rho _{\infty }}, \quad v_1 = -\frac {\textrm {i}P_1(x) \textrm {e}^{-\alpha _1 y_1}}{\rho _{\infty }}, \quad T_1 = 0 \quad \textrm {and} \quad \phi _1 = 0. \end{align}

Continuity of pressure requires

(5.9) \begin{align} P_1(x)=p_2(x) \quad \textrm {as} \quad y_1\rightarrow 0. \end{align}

Similarly, continuity of the wall-normal velocity $\tilde {v}$ between the main deck solution (5.4b ) and the upper deck solution (5.8c ) yields the condition

(5.10) \begin{align} \alpha _1 A(x) = \frac {p_2(x)}{\rho _{\infty }}. \end{align}

5.3. The lower deck

Recall that the concentration layer has a characteristic thickness of $O({\textit{Re}}^{-1/2}{\textit{Sc}}^{-1/3})$ . By setting ${\textit{Sc}}^{-1/3}\sim {\textit{Re}}^{-1/8}$ , the lower deck coincides with the concentration layer.

To match with the main deck, in the lower deck $y=\varepsilon ^5 y_3$ , for $y_3=O(1)$ . Perturbations in the lower deck are then expanded as

(5.11a–e) \begin{align} \tilde {u} &{}={} \left (u_3+O(\varepsilon )\right )E, \qquad \quad \tilde {v} {}={} \left (\varepsilon ^2 v_3+O\big(\varepsilon ^3\big)\right )E, \nonumber \\ \tilde {p} &{}={} \left (\varepsilon p_3+O\big(\varepsilon ^2\big)\right )E, \quad\,\, \tilde {T} {}={} \left (T_3+O(\varepsilon )\right )E, \nonumber \\ \tilde {\phi } &{}={} \left (\phi _3+O(\varepsilon )\right )E, \end{align}

where $u_3=u_3(x,y_3)$ , etc. As before, similar expansions are introduced for the perturbed quantities $\tilde {\mu }$ , $\tilde {\rho }$ , $\tilde {c}$ and $\tilde {k}$ .

In the main deck, the base velocity behaves as $U_{\!B}\sim \lambda y_2$ as $y_2\rightarrow 0$ , where $\lambda = U_{\!B,y_2}|_{y_2=0} (\equiv \rho _wf''(0)/x^{1/2})$ , and consequently from (5.4a ) and (5.4b )

(5.12a,b) \begin{align} u_2 \rightarrow \lambda A(x) \quad \textrm {and} \quad v_2 \rightarrow -\textrm {i}\alpha _1 \lambda A(x)y_2 \quad \textrm {as} \quad y_2 \rightarrow 0. \end{align}

Therefore, within the lower deck, the base flow is given by

(5.13a–d) \begin{align} U_{\!B} &{}={} \varepsilon \lambda y_3 + O\big(\varepsilon ^2\big), \quad &V_{\!B} &{}={} -\frac {1}{2}\varepsilon ^2\lambda _x y_3^2 + O\big(\varepsilon ^3\big), \nonumber \\ T_{\!B} &{}={} T_{w} + \varepsilon \sigma y_3 + O\big(\varepsilon ^2\big), \quad &\phi _B &{}={} \phi _{\infty } + \varepsilon \psi (x,y_3) + O\big(\varepsilon ^2\big), \end{align}

where $\sigma = T_{B,y_2}|_{y_2=0} (\equiv \rho _w\theta '(0)/x^{1/2})$ .

Substituting (5.11) and (5.13) into the linear stability equations (4.3) gives

(5.14) \begin{align} p_3=p_2(x), \end{align}

to match with the pressure in the main deck, and

(5.15a) \begin{align} u_3&=B(x)\int _{\chi _0}^\chi \textrm {Ai}(\grave {\chi }) \;\textrm {d}\grave {\chi }, \end{align}

(5.15b) \begin{align} p_2 &= -\frac {\omega _1\rho _{\infty }}{\alpha _1}\frac {B(x)\textrm {Ai}'(\chi _0)}{\chi _0}, \end{align}

where $B$ is an unknown, amplitude function, $\textrm {Ai}$ is the Airy function and

(5.16) \begin{align} \chi = \left (\frac {\textrm {i}\alpha _1\lambda \rho _{\infty }}{\mu _{\infty }}\right )^{1/3}\left (y_3 - \frac {\omega _1}{\alpha _1\lambda }\right )\!, \end{align}

for $\chi _0=\chi |_{y_3=0}$ .

Matching the streamwise velocity $\tilde {u}$ between the main deck solution (5.12a ) and the lower deck solution (5.15a ), gives

(5.17) \begin{align} B(x)\int _{\chi _0}^\infty \textrm {Ai}(\chi ) \; \textrm {d} \chi = \lambda A(x). \end{align}

Eliminating $A$ , $B$ and $p_2$ from (5.10), (5.15b ) and (5.17) yields the leading-order eigenrelation

(5.18a,b) \begin{align} \frac {\textrm {Ai}'(\chi _0)}{\int _{\chi _0}^\infty \textrm {Ai}(\chi ) \; \textrm {d}\chi } = \left (\frac {\textrm {i}\alpha _1\lambda \rho _{\infty }}{\mu _{\infty }}\right )^{1/3}\frac {\alpha _1}{\lambda ^{2}}, \end{align}

which, following the parameter scaling

(5.19a) \begin{align} \alpha _1 = \lambda ^{5/4}\left (\frac {\mu _{\infty }}{\rho _{\infty }}\right )^{1/4}\overline {\alpha } \quad \textrm {and} \quad \omega _1 = \lambda ^{3/2}\left (\frac {\mu _{\infty }}{\rho _{\infty }}\right )^{1/2}\overline {\omega }, \end{align}

becomes

(5.20a,b) \begin{align} \frac {\textrm {Ai}'(\chi _0)} {\int _{\chi _0}^\infty \textrm {Ai}(\chi ) \; \textrm {d}\chi } =\textrm {i}^{1/3}\overline {\alpha }^{4/3} \quad \textrm {for} \quad \chi _0=-\textrm {i}^{1/3}\frac {\overline {\omega }} {\overline {\alpha }^{2/3}}. \end{align}

For neutral stability, $\alpha _1, \alpha _2$ , etc. must be real, requiring $\chi _0 = -2.298 \textrm {i}^{1/3}$ and

(5.21) \begin{align} \frac {\textrm {Ai}'(\chi _0)} {\int _{\chi _0}^\infty \textrm {Ai}(\chi ) \; \textrm {d}\chi } =1.001\textrm {i}^{1/3}. \end{align}

Consequently, the neutral values of $\alpha _1$ and $\omega _1$ are given as

(5.22a) \begin{align} \alpha _1 & = 1.001\hat {\lambda }^{5/4}\left (\frac {\mu _{\infty }}{\rho _{\infty }}\right )^{1/4}x^{-5/8}, \end{align}

(5.22b) \begin{align} \omega _1 & = 2.299\hat {\lambda }^{3/2}\left (\frac {\mu _{\infty }}{\rho _{\infty }}\right )^{1/2}x^{-3/4}, \end{align}

where $\hat {\lambda }=\rho _wf''(0)$ . This gives the leading-order approximation for the frequency of the lower branch in terms of the Reynolds number $\textit{Re}$ :

(5.23) \begin{align} \omega _N \sim 2.299[\delta _1\hat {\lambda }]^{3/2}\left (\frac {\mu _{\infty }}{\rho _{\infty }}\right )^{1/2}Re^{-1/2}. \end{align}

Notably, in the limit ${\textit{Sc}}\rightarrow \infty$ , $\delta _1\hat {\lambda }\approx 0.572$ across all nanoparticle materials and $\phi _{\infty }$ . Thus, $2.299[\delta _1\hat {\lambda }]^{3/2}\approx 0.994$ .

Figure 16 depicts the gradient of the frequency $\omega _N$ , defined as $\Delta \omega _N=0.994[\mu _{\infty }/\rho _{\infty }]^{1/2}$ , as a function of $\phi _{\infty }$ for all seven nanoparticle materials listed in table 1. A gradient $\Delta \omega _N\lt 0.994$ indicates a destabilising effect, while $\Delta \omega _N\gt 0.994$ corresponds to stabilising behaviour. The solutions are qualitatively similar to and consistent with the linear stability results shown in figure 14(a): less dense materials are stabilising and denser materials are destabilising.

Figure 16.

Gradient $\Delta \omega _N=0.994[\mu _{\infty }/\rho _{\infty }]^{1/2}$ of the lower branch (5.23) as a function of $\phi _{\infty }$ for different nanoparticle materials.