Linear stability of nanofluid boundary-layer flow over a flat plate

Abstract

The linear stability of nanofluid boundary-layer flow over a flat plate is investigated using a two-phase formulation that incorporates the Brinkman (1952 J. Chem. Phys., vol. 20, pp. 571–581) model for viscosity along with Brownian motion (BM) and thermophoresis (TP), building upon the earlier work of Buongiorno (2006 J. Heat Transfer, vol. 128, pp. 240–250). Solutions to the steady boundary-layer equations reveal a thin nanoparticle concentration layer near the plate surface, with a characteristic thickness of $O({\textit{Re}}^{-1/2}{\textit{Sc}}^{-1/3})$, for a Reynolds number ${\textit{Re}}$ and Schmidt number ${\textit{Sc}}$. When BM and TP are neglected, the governing equations reduce to the standard Blasius formulation for a single-phase fluid, and the nanoparticle concentration layer disappears, resulting in a uniform concentration across the boundary layer. Neutral stability curves and critical conditions for the onset of the Tollmien–Schlichting (TS) wave are computed for a range of nanoparticle materials and volume concentrations. Results indicate that while the effects of BM and TP are negligible, the impact of nanoparticle density is significant. Denser nanoparticles, such as silver and copper, destabilise the TS wave, whereas lighter nanoparticles, like aluminium and silicon, establish a small stabilising effect. Additionally, the viscosity model plays a crucial role, with alternative formulations leading to different stability behaviour. Finally, a high Reynolds number asymptotic analysis is undertaken for the lower branch of the neutral stability curve.

1. Introduction

This paper is concerned with the influence of nanofluids on the linear stability of disturbances in the boundary-layer flow over a flat plate. Nanofluids are fluids containing nanoscale particles ranging from 1 to 100 nm, dispersed in a base fluid like water. These nanoparticles, composed of metal-based or carbon-based materials, enhance the thermal properties of the base fluid.

Since the seminal work of Choi (1995), nanofluids have received considerable interest, with a rapid growth in annual publications (Taylor et al. 2013). Numerous studies have investigated the thermal benefits of nanofluids, including comprehensive reviews by Das, Choi & Patel (2006), Wang & Mujumdar (2008a , b ), Kakaç & Pramuanjaroenkij (2009), Mahbubul, Saidur & Amalina (2012) and Mishra et al. (2014). These thermal improvements have led to a wide range of heat transfer applications, including cooling systems for automotive engines (Sidik, Yazid & Mamat 2015), electronics (Bahiraei & Heshmatian 2018), nuclear systems (Buongiorno & Hu 2009), solar thermal systems (Khullar et al. 2012), biomedical processes (Sheikhpour et al. 2020) and industrial applications (Wong & Leon 2010).

Despite the ongoing interest in nanofluids for their thermal benefits, relatively few investigations have examined the impact of nanofluids on the hydrodynamic stability of flows. This study aims to address this knowledge gap by investigating the capabilities of nanofluids in controlling laminar–turbulent transition processes.

1.1. Modelling nanofluid flows

A key aspect of modelling nanofluid flows is how suspended nanoparticles modify the fluid’s effective viscosity. For dilute suspensions of rigid, spherical particles, Einstein (1906) showed that the dynamic viscosity increases linearly with the nanoparticle volume concentration $\phi$ . He defined the effective dynamic viscosity as

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

where $\mu _{\textit{bf}}^*$ is the dynamic viscosity of the base fluid. Since Einstein’s work, many viscosity models have been proposed to account for additional factors, including particle shape, size distribution and particle–particle interactions. Batchelor (1977) extended Einstein’s formula to include the effects of Brownian motion (BM) (that is, the random movement of nanoparticles in a base fluid), while Brinkman (1952) proposed a semiempirical correlation valid for nanoparticle volume concentrations up to approximately $4\,\%$ . (The formulae for the Batchelor and Brinkman models are given in the subsequent section.) Comprehensive reviews of nanofluid viscosity models, including experimental and theoretical developments, are provided by Wang & Mujumdar (2008a ) and Mishra et al. (2014).

Another key aspect of nanofluid modelling is the treatment of the fluid either as a single-phase or a two-phase flow. Single-phase models treat the nanofluid as a homogeneous mixture with effective properties, while two-phase models account for interactions between the base fluid and nanoparticles. The latter approach can capture additional effects such as particle migration, BM and thermophoresis (TP) (that is, the movement of nanoparticles in a base fluid due to a temperature gradient). Moreover, two-phase flow models include a continuity equation for the nanoparticle volume concentration.

The steady boundary-layer flow over a flat plate has been investigated by Buongiorno (2006), Avramenko, Blinov & Shevchuk (2011) and MacDevette, Myers & Wetton (2014). These studies employed the Brinkman (1952) model to describe the nanofluid viscosity and incorporated BM and TP into the governing equations. To simplify the analysis, Buongiorno (2006) assumed the flow to be incompressible, even though modelling the nanofluid as a two-component mixture implies a non-constant density. Despite this apparent inconsistency, Buongiorno (2006) showed that the effects of BM and TP are negligible in nanofluids and attributed the observed heat transfer benefits to the improved thermophysical properties of the nanoparticles.

While acknowledging that BM and TP effects are weak, Avramenko et al. (2011) derived boundary-layer equations similar to those of Blasius (1908). However, despite accounting for compressibility effects in the base flow, the study implemented several simplifying assumptions. Notably, the incompressible flow condition was applied to the nanoparticle continuity equation (see equations (1)–(4) of Avramenko et al. (2011)). Additionally, the coefficients for BM and TP, defined below in (2.7), were treated as constants, even though they depend on temperature and nanoparticle volume concentration, respectively. Yet despite these simplifications, Avramenko et al. (2011) demonstrated that a thin concentration layer forms near the plate surface. This concentration layer modifies the velocity and temperature fields in the near-wall region, which may, in turn, influence instabilities within the boundary layer.

Both Buongiorno (2006) and Avramenko et al. (2011) confirmed that heat transfer, measured by the Nusselt number $Nu$ , is enhanced as the nanoparticle volume concentration $\phi$ increases. In contrast, MacDevette et al. (2014), who also confirmed that BM and TP are negligible, observed a reduced heat transfer coefficient as $\phi$ increases. They attributed the discrepancy with earlier studies to differences in the definition of the heat transfer coefficient.

The study of nanofluids in boundary-layer flows has been extended to include flows past vertical plates (Kuznetsov & Nield 2010), planar wall jets (Turkyilmazoglu 2016), stretching sheets (Reddy et al. 2025) and the flow due to a rotating-disk (Bachok, Ishak & Pop 2011; Turkyilmazoglu 2014; Mehmood & Usman 2018), with these studies reporting enhanced heat transfer due to the introduction of nanoparticles.

Using triple-deck theory, Wasaif (2023) modelled a nanofluid boundary-layer flow past a hump, on an otherwise flat plate. The study demonstrated that a nanofluid can suppress the region of flow separation along the rear side of the bump. More recently, Gandhi, Nepomnyashchy & Oron (2025) examined thermosolutal instabilities in a nanofluid layer with a deformable surface, showing how the Soret effect and thermal properties influence instability characteristics.

1.2. Linear stability studies

The linear stability of the incompressible Blasius boundary layer has been extensively studied, beginning with the seminal investigations of Tollmien (1933) and Schlichting (1933), which led to the Orr–Sommerfeld equation. These studies employed the parallel flow approximation, where the flow is assumed to be unidirectional and depends only on the wall-normal coordinate. The theoretical predictions for the Tollmien–Schlichting (TS) wave were subsequently confirmed experimentally by Schubauer & Skramstad (1947). Further theoretical and experimental insights into the stability of TS waves were reported by Jordinson (1970), Barry & Ross (1970), Ross et al. (1970) and Gaster (1974), amongst many others.

Using triple-deck theory, Smith (1979) undertook an asymptotic, high Reynolds number ${\textit{Re}}$ analysis to describe the structure of the lower branch of the neutral stability curve in the Blasius boundary layer (Lin 1955). The triple-deck framework consists of three layers: an upper deck, representing the inviscid outer flow and spans a thickness of $O({\textit{Re}}^{-3/8})$ ; a main deck, corresponding to the boundary layer, with thickness $O({\textit{Re}}^{-4/8})$ ; and a lower deck, a thin viscous sublayer of thickness $O({\textit{Re}}^{-5/8})$ , where viscous–inviscid interactions are dominant. (A formal definition for the Reynolds number ${\textit{Re}}$ is given below in (2.15a )) A subsequent study by Bodonyi & Smith (1981) employed a multideck approach to derive the corresponding structure of the upper branch of the neutral stability curve. Later, Smith (1989) extended the asymptotic analysis of the lower branch to compressible boundary-layer flows.

Building on earlier studies, Bertolotti, Herbert & Spalart (1992) employed parabolised stability equations to investigate both the linear and nonlinear development of TS waves in the Blasius boundary layer. Healey (1995) compared the asymptotic scalings of the lower and upper branches with solutions from the Orr–Sommerfeld equation and experimental observations. More recently, both asymptotic and numerical approaches have been utilised to model the effects of non-Newtonian viscosity (Griffiths, Gallacher & Stephen 2016) and temperature-dependent viscosity (Miller et al. 2018) on the stability of the Blasius boundary layer.

To the authors’ knowledge, there are only two previous studies concerning the linear stability of nanofluid boundary-layer flows. The first, by Turkyilmazoglu (2020), considered the application of nanofluids to several configurations, including the Kelvin–Helmholtz instability, Rayleigh–Bénard convection, instabilities in rotating disk flows and instabilities in the boundary-layer flow over a flat plate. Turkyilmazoglu modelled the latter flow as a single-phase flow, with quantities scaled on nanofluid properties, i.e. the combined characteristics of the base fluid and nanoparticles. This approach led to a Reynolds number based on nanofluid characteristics and a base flow described by the Blasius equation. The findings suggest that the Reynolds number of the nanofluid can be predicted using the Reynolds number of the base fluid. Moreover, results indicate that denser nanoparticle materials, like silver (Ag), stabilise the flow, while less dense nanoparticle materials, such as alumina ( $\text{Al}_2\text{O}_3$ ), destabilise the flow at sufficiently larger volume concentrations $\phi$ . However, the rationale for scaling quantities on nanofluid characteristics is unclear, as the resulting Reynolds number changes as the nanoparticle volume concentration $\phi$ increases, making it difficult to compare solutions. In the following study, the nanofluid flow is modelled as a two-phase flow that includes diffusion effects due to BM and TP, with the Reynolds number based on the base fluid properties to facilitate comparisons across different nanoparticle materials and variable $\phi$ .

A second study, by Laouer et al. (2024), examined the linear stability of a nanofluid flow past stationary and moving wedges. Similar to Turkyilmazoglu (2020), Laouer et al. (2024) employed a single-phase flow approach, using the base flow formulation of Yacob, Ishak & Pop (2011) and a linear stability analysis that simplifies to the standard Orr–Sommerfeld equation for a regular fluid. Laouer et al. (2024) showed that, for a nanofluid flow over a stationary wedge due to a favourable pressure gradient, increasing the volume concentration $\phi$ leads to a destabilising effect. Additionally, Laouer and coworkers suggest that heavier nanoparticle materials, such as copper (Cu), have a stabilising effect, while lighter materials, like titanium oxide ( $\text{TiO}_2$ ) and $\text{Al}_2\text{O}_3$ , destabilise the flow. However, this latter finding appears to contradict the results presented in figure 8 of their paper, which shows that Cu nanoparticles shift neutral stability curves to the left and smaller Reynolds numbers, while $\text{TiO}_2$ and $\text{Al}_2\text{O}_3$ nanoparticles shift neutral stability curves to the right and higher Reynolds numbers.

1.3. Outline of paper

The following study investigates the linear stability of nanofluid flow over a flat plate using a two-phase flow model that accounts for BM and TP. This model addresses the inconsistencies in previous single-phase studies and provides a physically consistent method for analysing stability trends. Both numerical and asymptotic analyses are undertaken to compute neutral stability curves and examine the lower branch behaviour at high Reynolds numbers. The most amplified TS disturbances appear near the lower branch of the neutral curve, and this, combined with the need to validate our numerical solutions, motivates the analysis of the lower rather than the upper branch.

The remainder of this paper is outlined as follows. The governing equations are introduced in the next section, followed by the steady, two-dimensional boundary-layer equations and its solutions in § 3. Linear stability results for three-dimensional disturbances, including neutral stability curves and critical conditions, are presented in § 4. An asymptotic analysis of the lower branch is provided in § 5. Conclusions are given in § 6.

2. Governing equations

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 2014), along with a continuity equation for the nanoparticles (Buongiorno 2006; Avramenko et al. 2011; MacDevette et al. 2014),

(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 (1881) 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 (2008a ).

The dynamic viscosity of the nanofluid $\mu ^*$ , used throughout the subsequent study, is given by the Brinkman (1952) 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. 2014). 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 (2008a ) and Mishra et al. (2014), which encompass properties such as the size and distribution of nanoparticles. For instance, Batchelor (1977) 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 (1998) and Maiga et al. (2004) 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 (1973), Buongiorno (2006) and MacDevette et al. (2014).

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. (2011), 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 (1952) 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 (1977), Pak & Cho (1998) and Maiga et al. (2004) models.

Table 1. Thermophysical properties of water and various materials used for nanoparticles, as given in Buongiorno (2006), Wang & Mujumdar (2008a ), Bachok et al. (2011), MacDevette et al. (2014), Turkyilmazoglu (2014, 2020) 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 (1952), Batchelor (1977), Pak & Cho (1998) and Maiga et al. (2004) 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. Steady boundary-layer flow

3.1. Boundary-layer equations

Following the derivation of Ruban (2017), 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. (2011), 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 (1995) 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. Linear stability analysis

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 1984; Schmid & Henningson 2001). 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 (2000). 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 1976; Grosch & Salwen 1978; Salwen & Grosch 1981; Schmid & Henningson 2001), the eigenspectrum consists of multiple branches. A discrete set of modes are located on the A-branch (Mack 1976) 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 2001).

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 (1952) 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. Asymptotic analysis

To describe the lower-branch structure of the neutral stability curve, we follow the approach of Smith (1979) and assume a large Reynolds number ${\textit{Re}}$ . Consequently, linear disturbances on the lower branch are governed by a triple deck structure with a main deck of thickness $O({\textit{Re}}^{-1/2})$ , an upper deck of thickness $O({\textit{Re}}^{-3/8})$ and a lower deck of thickness $O({\textit{Re}}^{-5/8})$ , with streamwise length $O({\textit{Re}}^{-3/8})$ and frequency $O({\textit{Re}}^{-1/4})$ . A diagram of the triple deck structure is shown in figure 15 for $\varepsilon ={\textit{Re}}^{-1/8}$ . In addition,

(5.1a,b) \begin{align} x=1+\varepsilon ^3 X \quad \textrm {and} \quad t=\varepsilon ^2\hat {t}, \end{align}

while linear disturbances are taken to be proportional to

(5.2a) \begin{align} E=\exp \left (\textrm {i}\left (\varTheta (X) - \omega \hat {t}\right )\right )\!, \end{align}

for

(5.2b,c) \begin{align} \frac {{\rm d}\theta }{{\rm d}X} = \alpha _1(x) +\varepsilon \alpha _2(x) + {\cdots} \quad \textrm {and} \quad \omega = \omega _1+\varepsilon \omega _2+{\cdots} . \end{align}

Figure 15. Diagram of the triple deck structure of the lower-branch of the neutral stability curve for $\varepsilon ={\textit{Re}}^{-1/8}$ . Regions $1$ , $2$ and $3$ correspond to the upper, main and lower decks, respectively.

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.

6. Conclusions

A linear stability study has been conducted on the nanofluid boundary-layer flow over a flat plate, extending the earlier work of Buongiorno (2006), Avramenko et al. (2011), MacDevette et al. (2014) and Turkyilmazoglu (2020). The model employs a two-phase flow formulation that incorporates the effects of BM and TP, with all quantities scaled on the base fluid characteristics, providing a physically consistent approach for investigating stability trends. Although the influence of BM and TP is relatively weak, a thin concentration layer with a characteristic thickness of $O({\textit{Re}}^{-1/2}{\textit{Sc}}^{-1/3})$ develops within the boundary layer, which modifies the near-wall velocity and temperature fields. The concentration layer disappears when BM and TP are ignored, with the nanoparticle volume concentration $\phi$ uniform throughout the boundary layer.

In terms of thermodynamic performance, all seven materials modelled herein establish an increasing Nusselt number $Nu$ , with greater benefits obtained for denser materials like Ag and Cu.

Despite the emergence of a thin concentration layer, numerical and asymptotic stability calculations show that BM and TP have a negligible impact on the onset of the TS wave. In fact, linear stability characteristics and the onset of TS waves can be accurately predicted using solutions to the standard Blasius flow, which effectively models the nanofluid flow as a single-phase fluid. The Reynolds number for the nanofluid is given as

(6.1) \begin{align} {\textit{Re}} = \frac {\mu \widehat {\textit{Re}}}{\rho }, \end{align}

for the Blasius flow Reynolds number $\widehat {\textit{Re}}$ . Consequently, the stability of the nanofluid boundary-layer flow is governed by the density $\rho$ and viscosity $\mu$ of the nanofluid. In particular, the density ratio $\hat {\rho }=\rho ^*_{\textit{np}}/\rho ^*_{\textit{bf}}$ is critical to determining whether the nanofluid is stabilising or destabilising. Denser nanoparticle materials, such as Ag and Cu, significantly destabilise the TS wave. In contrast, a small stabilising effect is achieved by lighter materials, like Al and Si. This observation differs from the one-phase flow study conducted by Turkyilmazoglu (2020), which predicted the opposite outcome. However, in Turkyilmazoglu’s investigation, physical quantities were scaled on the characteristics of the nanofluid rather than the base fluid, leading to a Reynolds number that varied with the type of nanoparticle material and volume concentration.

The results presented above are based on a nanofluid with water as the base fluid. Replacing water with a less dense fluid, like ethanol, would increase the density ratio $\hat {\rho }$ for all materials. While this change would enhance the thermal benefits of the nanofluid, it would lead to a further destabilisation of the TS wave, even for those nanofluids composed of lighter materials like Al and Si.

Another key factor influencing nanofluid stability is the choice of viscosity model. The above study adopted the Brinkman (1952) model (2.5) to represent the nanofluids dynamic viscosity, ensuring consistency with earlier investigations. However, alternative models can produce very different results. For instance, the correlations due to Pak & Cho (1998) and Maiga et al. (2004) in pipe and tube flows (see (2.6b )) predict larger increases in viscosity as the nanoparticle volume concentration $\phi$ increases. Assuming these models can be applied to the boundary-layer flow on a flat plate, the corresponding stability calculations based on (6.1) indicate a strong stabilising effect for all nanoparticle materials, in contrast to the destabilising trends observed for the Brinkman model. Since the Reynolds number and resulting stability characteristics depend on the viscosity model, accurately determining $\mu$ is crucial. The variety of models summarised in Wang & Mujumdar (2008a ) and Mishra et al. (2014), including those that include nanoparticle aggregation, size and shape effects, and temperature-dependent viscosity, highlights the need for further experimental measurements of nanofluid viscosity in boundary-layer flows to enable the selection of the correct model and ensure physically accurate stability predictions.

In addition to a base fluid of water, nanoparticles were assumed to have a diameter of $d_{\textit{np}} = 20$ nm with a free stream temperature $T_{\infty } = 300$ K. Varying the nanoparticle diameter over the range $d_{\textit{np}} = 1$ nm to $d_{\textit{np}} = 100$ nm yields Schmidt numbers ${\textit{Sc}}$ ranging from approximately $O(10^3)$ to $O(10^5)$ . Even at the lower end of this range, ${\textit{Sc}}$ is sufficiently large that BM and TP effects remain negligible. Furthermore, applying the analysis to an ethanol-based fluid establishes comparable ${\textit{Sc}}$ values, whereas an oil-based fluid results in even larger values due to its higher viscosity. While increasing $T_{\infty }$ can reduce ${\textit{Sc}}$ to $O(10^2)$ for $T_{\infty } = 1000$ K, such high temperatures are generally unrealistic for practical nanofluid applications.

Future investigations into nanofluid boundary-layer flows could include non-parallel and nonlinear stability effects by using parabolised stability equations, following the approach of Bertolotti et al. (1992). Additionally, the analysis may be applied to other geometries, including rotating disk boundary layers and wall jets, as considered by Turkyilmazoglu (2020). However, based on the observations above, we anticipate that the stability trends would remain similar: heavier nanoparticles are expected to be destabilising, lighter nanoparticles stabilising, with the flows well-approximated by the standard Blasius flow without nanoparticles.