3.3.1. NSC and mode of instability: role of spanwise wavenumber and slip length

The impact of spanwise wavenumber on neutral stability curves in the $(\alpha , \textit{Re})$ -plane for $\hat{d}=0.1$ is shown in figure 4(a–c) for $L_x=0$ and figure 4(d–f) for $L_x=0.02$ . Six values of $\beta$ , including two critical values 0 and 1.35 for $L_x$ equal to 0 and 0.02, respectively, are chosen. The corresponding wave propagation speed in the $(\alpha , c_r)$ -plane is illustrated by the dashed lines that assist in identifying the instability mode. The CIP for $L_x$ equal to 0 and 0.02 are $(6322.50, \, 0.808, \, 0)$ and $(5188.39, \, 1.496, \, 1.35)$ , respectively. In figures 4(a) and 4(d), each NSC shows two minima in two different lobes. The minimum in the left lobe appears at relatively small (large) values of $\alpha$ ( $c_r$ ) and corresponds to the porous mode of instability (Chang et al. Reference Chang, Chen and Straughan2006). The second minimum appears in the right lobe at relatively large (small) values of $\alpha$ ( $c_r$ ) and corresponds to the fluid mode of instability. In each NSC, the 2-D critical point is marked by a singleton bullet. The profile of wave speed shows a monotonic decrease with $\alpha$ in the first lobe and remains almost constant in the second lobe except in the vicinity of the end part of the second lobe, where $\textit{Re}$ increases rapidly and $c_r$ decreases on increasing the value of $\alpha$ . As $\beta$ is increased to $0.5$ , the instability mode at the 2-D critical point shifts from porous to fluid on changing the value $0$ of $L_x$ by $0.02$ . For $\beta = 1$ , the NSC at $L_x$ equal to $0 \, (0.02)$ possesses two (three) lobes, left and right (left, central and right), that correspond to two (three) different modes of instability, which will be analysed shortly. The minimum in the left (central) lobe dominates the instability.

In figure 4(e), the profile of wavespeed experiences a jump discontinuity as it moves from the left to the central lobe and shows a possibility of a sudden change of mode of instability. A slight increase in $\beta$ to $1.1$ results in an expansion of the unstable region in the central lobe, accompanied by a reduction in the unstable areas within the left and right lobes, and the 2-D critical point shifts downward. The left lobe vanishes completely and the minima in the central lobe become the least one for $\beta =1.35$ (see figure 4 f). Thus, for $L_x=0.02$ , the three-dimensional disturbances are least stable. Similar behaviour was also observed in the linear stability of compressible boundary layer flows with non-ideal fluids (Ren, Marxen & Pecnik Reference Ren, Marxen and Pecnik2019; Ren & Kloker Reference Ren and Kloker2022; Ren et al. Reference Ren, Wu, Mao, Wang and Kloker2025), where the appearance of three-dimensional perturbations as the least stable is possible in the transcritical and supercritical flow regimes. A further increase of $\beta$ to $1.5$ causes the 2-D critical point to rise again and reduces the right lobe’s regime, confirming the surpassing of the most unstable spanwise wavenumber. Thus, in contrast to $L_x=0$ , where the 2-D critical point lies in either the left or right lobe, the same for $L_x=0.02$ based on the value of $\beta$ may lie in the left, right or central lobe.

Figure 5.

Normalised disturbed velocity contour in the wall-normal direction at 2-D critical point, for $L_x=0$ when (a) $\alpha =0.808$ , $\beta =0$ , $\textit{Re}=6322.51$ ; (b) $\alpha = 2.063$ , $\beta =0.5$ , $\textit{Re}=8960.87$ ; (c) $\alpha =1.90$ , $\beta =1.35$ , $\textit{Re}=9951.38$ ; and for $L_x=0.02$ when (d) $\alpha =0.767$ , $\beta =0$ , $\textit{Re}=6559.38$ ; (e) $\alpha =2.241$ , $\beta =0.5$ , $\textit{Re}=7650.00$ ; ( f) $\alpha =1.496$ , $\beta =1.35$ , $\textit{Re}=5188.38$ .

The disturbance velocity contour plots are used to explicitly understand the mode of instability in different lobes. For this, three 2-D critical points for $L_x=0$ ( $L_x=0.02$ ) corresponding to $\beta =0$ , $0.5$ and $1.35$ are chosen. The flow patterns for the normalised disturbed velocities ( $u_{\!p}$ and $u$ ) in the $xz$ -plane for one wave period (i.e. $x \in [0,2\pi /\alpha ]$ ) are displayed in figure 5(a–c) (figure 5 d–f) at the respective above-mentioned three points. As shown in figures 5(a) and 5(d), when the instability dominates in the left lobe, the disturbance initiates in the fluid region and penetrates through the interface into the porous region. This results in an advective flow inside the porous layer, which is mainly responsible for the onset of instability. In this case, the instability mode is defined as the porous mode (PM) of instability. In figures 5(b), 5(c), 5(e) and 5( f), the flow patterns are mostly confined to the fluid region. In figures 5(b), 5(c) and 5(e), when the instability lies in the right lobe, the contours are tilted, and there is a limited penetration of disturbance inside the porous medium. However, when the instability lies in the central lobe, the contours exhibit square-like patterns in figure 5( f) and have no penetration inside the porous medium. The pattern of eigenfunctions corresponding to figures 5(b), 5(c), 5(e) and 5( f) reveals the anti-symmetry and symmetry about $z=0.5$ , respectively. Thus, following Chang et al. (Reference Chang, Chen and Straughan2006, Reference Chang, Chen and Chang2017) and Sharma et al. (Reference Sharma, Aleria and Bera2024), the modes of instabilities in the former and later cases are defined as odd-fluid mode (OFM) and even-fluid mode (EFM), respectively. The mode of instability in the first and second lobes in figures 4(a), 4(b), 4(d) (figure 4 f) is also checked and found as PM and OFM (EFM and OFM), respectively. However, the same in the left, central and right lobes in figure 4(e) are identified as PM, EFM and OFM, respectively. It has also been checked that a jump in wave speed profile in figure 4(e) leads to a sudden change of mode of instability from porous to even-fluid.

Figure 6.

Neutral stability curves for ${\textit{Da}}=10^{-6}$ at $\hat{d}=0.1$ : (a) $\beta =0$ , (b) $\beta =0.98$ . The dashed curves represent the corresponding wave propagation speeds. The insets in panels (a) and (b) represent the normalised disturbed velocity profile at $(\alpha , \beta , \textit{Re})= (0.65, 0, 7311.84)$ (dashed line), $(\alpha , \beta , \textit{Re})= (1.972, \,0, \, 3334.90)$ (solid line) and $(\alpha , \beta , \textit{Re})= (1.889, \,0.98, \, 3950.86)$ (solid line) for $L_x=0.05$ .

So, from the above analysis, we can see that for non-zero slip length, depending on the value of $\beta$ , the mode of instability along the NSC may be of PM, OFM or EFM. The PM and OFM vanish for relatively higher values of $\beta$ . To check whether the mode of instability for $\beta$ equal to 0 or close to 1 can be free from OFM for some higher values of $L_x$ or not, the NSC is also plotted for two different values 0.05 and 0.08 of $L_x$ at 0 and 0.98 of $\beta$ in figures 6(a) and 6(b), respectively. Note that the CIP for $L_x$ equal to 0.05 and 0.08 are $(1.889, \, 0.98, \, 3950.86)$ and $(1.972, \, 0, \, 3334.90)$ , respectively. From figure 6(a), it can be seen that the NSC exhibits bimodal characteristics and the wave speed experiences a sudden change upon shifting from one lobe to another. Here, it has also been checked that along the NSC, the mode of instability changes from porous to even-fluid (see the inset in figure 6 a). However, the inset in figure 6(b) shows an even symmetry of the normalised disturbed velocity profile concerning $z=0.5$ . Thus, the unimodal behaviour of NSCs for $\beta =0.98$ demonstrates the even-fluid mode of instability.

3.3.2. Effect of slip length on 2-D critical point as a function of $\beta$ for different values of $\hat{d}$

Thus far, the instability characteristics along the NSCs in figures 4 and 6 were investigated for a fixed depth ratio $\hat{d}=0.1$ and four discrete values of $L_x$ . For $L_x$ equal to zero, the instability mode at CIP is porous; for the other three values, it is even-fluid. The critical value of $\textit{Re}$ at $L_x=0$ is higher than that at $L_x = 0.08$ . Considering four values (0, 0.02, 0.05, 0.08) of $L_x$ , we can see that the least stable disturbances change from two-dimensional to three-dimensional and then again to two-dimensional. The change of dimension of the least stable disturbances on changing the value of slip length makes us curious to know that, under such slip conditions, how do the same changes on varying the values of $\hat{d}$ ? What will be the role of $\beta$ on it? To answer these queries and shed more light on the mechanism of instability at each 2-D critical point, figures 7(a–c) and 7(d–f) are plotted for $L_x$ equal to $0.02$ and $0.05$ , respectively. The dot point in each curve represents the location of $\beta _c$ .

Figure 7.

Variation of 2-D critical parameters $\textit{Re}$ , $\alpha$ and $c_r$ with $\beta$ for different values of $\hat{d}$ , when (a,b,c) $L_x=0.02$ , and (d,e, f) $L_x=0.05$ . The dot represents the location of $\beta _c$ for the critical instability point.

In figure 7(a–c), when $\hat{d}$ is equal to 0.08 and 0.5, the two-dimensional disturbances dominate the flow instability, i.e. CIP is two-dimensional; however, for other values of $\hat{d}$ , it is three-dimensional. To identify the instability mode at any 2-D critical point, we shall follow the definition of the mode of instability as a function of $c_r$ and by looking at the nature of the corresponding eigenfunction. As discussed in the previous subsection along the NSC, for a porous (even-fluid) mode of instability, the value of $c_r$ is high (low) and the $c_r$ profile experiences a jump whenever the PM instability is changed into EFM. When the value of $c_r$ continuously decreases and shows a constant nature for a small sub-interval of $\alpha$ and then decreases further, it indicates the change of instability from PM to OFM. However, when EFM is changed into OFM, it introduces either a jump or a cusp point on the profile of $c_r$ . It has been checked that a similar characteristic of the $c_r$ profile as a function of $\beta$ can be correlated for identifying the mode of instability for any 2-D critical point.

Due to the relatively small (large) value of $\alpha _c$ as well as the large (small) value of $c_r$ at $\hat{d}$ equal to 0.08 (0.5), the mode of instability at the critical instability point (CIP) is porous (even-fluid), which is also confirmed by checking the nature of corresponding eigen-function in figure 21 (see Appendix B). For $L_x = 0.02$ , we have already seen that the mode of instability at CIP is EFM for $\hat{d}=0.1$ , which is also the mode of instability at CIP for other depth ratios, except for $\hat{d}=0.08$ (see figure 21). In contrast, for $L_x=0.05$ , the mode of instability at CIP for all considered values of $\hat{d}$ is EFM. Furthermore, the instability boundary curve at $L_x = 0.02$ , $\hat{d}=0.08$ have two lobes. In the first lobe up to $\beta =1.45$ , the value $\alpha$ ( $c_r$ ) increases (decreases) smoothly, and then jumps to a large value and remains almost constant (jumps to a small value and remains practically constant). The second lobe starts with a jump discontinuity, where the value of $\alpha$ ( $c_r$ ) increases near the discontinuity and then decreases smoothly (experiences a sudden drop and then increases). It has been checked that in the first lobe, the mode of instability is PM except in the sub-interval of $\beta$ in which $\alpha$ experiences a jump and remains almost constant, and the mode of instability is checked as OFM. In the second lobe, $\alpha$ ( $c_r$ ) decreases (increases) smoothly and the mode of instability of any 2-D critical point is found as EFM. For $\hat{d}=0.1$ , the instability modes at $\beta$ equal to 0 and 0.5 are checked as PM and OFM, respectively, whereas at $\beta$ equal to $1, \, 1.1, \, 1.35$ and $1.5$ are EFM (see figure 4 d–f). For other values of $\hat{d}$ , EFM prevails as the mode of instability along the entire 2-D instability boundary.

As the value, $0.02$ of $L_x$ is replaced by 0.05, only at $\hat{d} = 0.08$ , the 2-D instability curve in figure 7(d) has two lobes. The first lobe shows a tendency to die out when the value of $L_x$ is replaced by 0.05. The profiles of wavenumber as well as wave speed for $L_x$ equal to 0.05 show the mode of instability at any 2-D critical point in the first lobe as PM, and in the second lobe as EFM.

Figure 8.

Effect of $L_x$ on (a) $\textit{Re}_c$ , (b) $\beta _c$ , (c) $\alpha _c$ and (d) $c_r$ for different values of $\hat{d}$ .

3.3.3. Three-dimensional critical point as a function of slip length, for different values of  $\hat{d}$

To understand the explicit role of anisotropic slip length on the CIP, the instability boundary curves in the ( $L_x, \textit{Re}_c$ )-, ( $L_x, \beta _c$ )-, ( $L_x, \alpha _c$ )-planes and the corresponding wave speed profile in the ( $L_x, c_r$ )-plane, with discretisation step $0.001$ for $L_x$ , are plotted in figure 8(a–d), respectively, for different values of $\hat{d}$ . It can be seen from figure 8(a–c) that each instability boundary curve can be split into three components. In the first component, i.e. in the vicinity of the no-slip condition, the least stable disturbance is two-dimensional, whereas in the second and third components, the least stable disturbances for each CIP are three-dimensional and two-dimensional, respectively. It can also be seen that the variation of $\textit{Re}_c$ with $L_x$ in the first component is primarily influenced by the critical mode of instability of CIP under the no-slip conditions. For $\hat{d} \lt 0.12$ , the PM instability dominates, while for $\hat{d} \gt 0.12$ , the EFM mainly triggers the instability under the no-slip case (see figure 8 d). Under the same conditions, for $\hat{d} = 0.12$ , the OFM emerges as the most unstable mode, which is confirmed by looking at the structure of the corresponding eigenfunction. In general, when PM or EFM dominates the instability under the no-slip condition, figure 8(a) indicates that a slight increase in $L_x$ from the no-slip case enhances the stability of the flow. In contrast, the same destabilises the flow when the OFM dominates the flow instability, as shown in the inset of figure 8(a). Based on the corresponding $c_r$ profile, we see that the instability mode in the second and third components is EFM. Furthermore, the length of the subinterval of $L_x$ in the second component decreases on increasing the depth of the fluid layer, i.e. by increasing the value of $\hat{d}$ , which is also expected since the least stable disturbances in the limiting-channel plane Poiseuille flow are two-dimensional.

4.1.1. Equations for different harmonics

After substituting the Fourier expansions (4.1)–(4.2) into the governing equations (2.1)–(2.2) and separating the harmonic components, the equations for different harmonics are given as follows:

$E^0$ -harmonic equations:

(4.5) \begin{align} \mathcal{L}\! \left ( U_0, P_0 \right ) + c_i | B |^2 \left ( \mathcal{L}\! \left ( U_1, P_1 \right ) - \textit{Re} \frac {\text{d}}{\text{d}z} \left (u_{10} \tilde {w}_{10}+\tilde {u}_{10} w_{10} \right )\right ) & = {O} \big(c_i^2 \big), \end{align}

(4.6) \begin{align} \mathcal{L}^p\! \left ( U_{0p}, P_{0p}, Da \right )+c_{ip} |B_{\!p}|^2 \mathcal{L}^p\! \left ( U_{1p}, P_{1p}, Da \right ) & = {O} \big(c_{ip}^{2} \big). \end{align}

$E^1$ -harmonic equations:

(4.7) \begin{align} & c_{i}^{{1}/{2}} B {\mathscr{L}}_{o}(\alpha , \beta , u_{10}, v_{10}, w_{10}) + c_{i}^{{3}/{2}} | B |^2 B {\mathscr{L}}_{o}(\alpha , \beta , u_{11}, v_{11}, w_{11})={O} \big (c_{i}^{{5}/{2}}\big ), \end{align}

(4.8) \begin{align} & c_{i}^{{1}/{2}} B {\mathscr{L}}_{x}(\alpha , \beta , c, p_{10}, u_{10}, w_{10}, U_{0}, \textit{Re})+ c_{i}^{{3}/{2}} | B |^2 B {\mathscr{L}}_{x}(\alpha , \beta , c, p_{11}, u_{11}, w_{11}, U_{0}, \textit{Re}) \nonumber \\ &\qquad +c_i^{{3}/{2}} \left [ \textit{Re} \frac {\text{d}B}{\text{d}\tau } u_{10} - \alpha B \textit{Re} u_{10} -B | B |^2 \mathcal{G}_x \right ]={O} \big (c_{i}^{{5}/{2}}\big ), \end{align}

(4.9) \begin{align} & c_{i}^{{1}/{2}} B {\mathscr{L}}_{y}(\alpha , \beta , p_{10}, v_{10}, U_{0}, \textit{Re})+ c_{i}^{{3}/{2}} | B |^2 B {\mathscr{L}}_{y}(\alpha , \beta , p_{11}, v_{11}, U_{0}, \textit{Re}) \nonumber \\ &\qquad +c_i^{{3}/{2}} \left [ \textit{Re} \frac {\text{d}B}{\text{d}\tau } v_{10} - \alpha B \textit{Re} v_{10} -B | B |^2 \mathcal{G}_y \right ]={O} \big (c_{i}^{{5}/{2}}\big ), \end{align}

(4.10) \begin{align} & c_{i}^{{1}/{2}} B {\mathscr{L}}_{z}(\alpha , \beta , p_{10}, w_{10}, U_{0}, \textit{Re})+ c_{i}^{{3}/{2}} | B |^2 B {\mathscr{L}}_{z}(\alpha , \beta , p_{11}, w_{11}, U_{0}, \textit{Re}) \nonumber \\ &\qquad +c_i^{{3}/{2}} \left [ \textit{Re} \frac {\text{d}B}{\text{d}\tau } w_{10} - \alpha B \textit{Re} w_{10} -B | B |^2 \mathcal{G}_z \right ]={O} \big (c_{i}^{{5}/{2}}\big ), \end{align}

(4.11) \begin{align} & c_{ip}^{{1}/{2}} B_{\!p} {\mathscr{L}}_{o}^{p}(\alpha _{\!p}, \beta _{\!p}, u_{10p}, v_{10p}, w_{10p})+ c_{ip}^{{3}/{2}} |B_{\!p}|^2 B_{\!p} {\mathscr{L}}_{o}^{p}(\alpha _{\!p}, \beta _{\!p}, u_{11p}, v_{11p}, w_{11p}) \nonumber \\ & \qquad \qquad ={{O}} \big (c_{ip}^{{5}/{2}}\big ), \end{align}

(4.12) \begin{align} & c_{ip}^{{1}/{2}} B_{\!p} {\mathscr{L}}_{x}^{p}(\alpha _{\!p}, c_{\!p}, p_{10p}, u_{10p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad+ c_{ip}^{{3}/{2}} |B_{\!p}|^2 B_{\!p} {\mathscr{L}}_{x}^{p}(\alpha _{\!p}, c_{\!p}, p_{11p}, u_{11p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad +\frac { \textit{Re}_{\!p} c_{ip}^{{3}/{2}}}{\epsilon } \left [\frac {\text{d}B_{\!p}}{\text{d}\tau _{\!p}} u_{10p} - \alpha _{\!p} B_{\!p} u_{10p} \right ]={{O}} \big (c_{ip}^{{5}/{2}}\big ), \end{align}

(4.13) \begin{align} & c_{ip}^{{1}/{2}} B_{\!p} {\mathscr{L}}_{y}^{p}(\alpha _{\!p}, \beta _{\!p}, c_{\!p}, p_{10p}, v_{10p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad + c_{ip}^{{3}/{2}} B_{\!p} |B_{\!p}|^2 {\mathscr{L}}_{y}^{p}(\alpha _{\!p}, \beta _{\!p}, c_{\!p}, p_{11p}, u_{11p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad +\frac { \textit{Re}_{\!p} c_{ip}^{{3}/{2}}}{\epsilon } \left [\frac {\text{d}B_{\!p}}{\text{d}\tau _{\!p}} v_{10p} - \alpha _{\!p} B_{\!p} v_{10p} \right ]={{O}} \big (c_{ip}^{{5}/{2}}\big ), \end{align}

(4.14) \begin{align} & c_{ip}^{{1}/{2}} B_{\!p} {\mathscr{L}}_{z}^{p}(\alpha _{\!p}, c_{\!p}, p_{10p}, u_{10p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad +\, c_{ip}^{{3}/{2}} |B_{\!p}|^2 B_{\!p} {\mathscr{L}}_{z}^{p}(\alpha _{\!p}, c_{\!p}, p_{11p}, u_{11p}, \textit{Re}_{\!p}, Da, \epsilon ) \nonumber \\ & \qquad \qquad +\, \frac { \textit{Re}_{\!p} c_{ip}^{{3}/{2}}}{\epsilon } \left [\frac {\text{d}B_{\!p}}{\text{d}\tau _{\!p}} w_{10p} - \alpha _{\!p} B_{\!p} w_{10p} \right ] ={{O}} \big (c_{ip}^{{5}/{2}}\big ). \end{align}

$E^2$ -harmonic equations:

(4.15) \begin{align} c_{i} B^2 & {\mathscr{L}}_{o} (2\alpha , 2\beta , u_{20}, v_{20}, w_{20}) ={{O}} \big (c_{i}^{2}\big ), \end{align}

(4.16) \begin{align} c_{i} B^2 & \Bigg [ {\mathscr{L}}_{x} (2\alpha , 2\beta , c, p_{20}, u_{20}, w_{20}, U_0, \textit{Re}) \nonumber \\ & + \textit{Re} \Bigg ( \text{i} \alpha u^{2}_{10} + \text{i} \beta u_{10} v_{10} + w_{10} \frac {\text{d} u_{10}}{\text{d} z} \Bigg ) \Bigg] ={{O}} \big (c_{i}^{2}\big ), \end{align}

(4.17) \begin{align} c_{i} B^2 & \Bigg [ {\mathscr{L}}_{y} (2\alpha , 2\beta , c, p_{20}, v_{20}, U_0, \textit{Re}) \nonumber \\ & + \textit{Re} \Bigg ( \text{i} \alpha u_{10} v_{10} + \text{i} \beta v^2_{10} + w_{10} \frac {\text{d} v_{10}}{\text{d} z} \Bigg ) \Bigg ] ={{O}} \big (c_{i}^{2}\big ), \end{align}

(4.18) \begin{align} c_{i} B^2 & \Bigg [ {\mathscr{L}}_{z} (2\alpha , 2\beta , c, p_{20}, w_{20}, U_0, \textit{Re}) \nonumber \\ & + \textit{Re} \Bigg ( \text{i} \alpha u_{10} w_{10} + \text{i} \beta v_{10} w_{10} + w_{10} \frac {\text{d} w_{10}}{\text{d} z} \Bigg ) \Bigg ] ={{O}} \big (c_{i}^{2}\big ), \end{align}

(4.19) \begin{align} c_{ip} B_{p}^2 & {\mathscr{L}}_{o}^{p} (2\alpha _{\!p}, 2\beta _{\!p}, u_{20p}, v_{20p}, w_{20p}) ={{O}} \big (c_{ip}^{2}\big ), \end{align}

(4.20) \begin{align} c_{ip} B_{p}^2 & {\mathscr{L}}_{x}^{p}(2\alpha _{\!p}, c_{\!p}, p_{20p}, u_{20p}, \textit{Re}_{\!p}, Da, \epsilon ) ={{O}} \big (c_{ip}^{2}\big ), \end{align}

(4.21) \begin{align} c_{ip} B_{p}^2 & {\mathscr{L}}_{y}^{p}(2\alpha _{\!p}, 2\beta _{\!p}, c_{\!p}, p_{20p}, v_{20p}, \textit{Re}_{\!p}, Da, \epsilon ) ={{O}} \big (c_{ip}^{2}\big ), \end{align}

(4.22) \begin{align} c_{ip} B_{p}^2 & {\mathscr{L}}_{z}^{p}(2\alpha _{\!p}, c_{\!p}, p_{20p}, w_{20p}, \textit{Re}_{\!p}, Da, \epsilon ) ={{O}} \big (c_{ip}^{2}\big ). \end{align}

The $\mathcal{G}$ functions defined in (4.8)–(4.10) are given by

(4.23) \begin{align} \mathcal{G}_x &=- \textit{Re} \biggl [ \text{i} \alpha \left ( u_{10} U_1 + \tilde {u}_{10} u_{20} \right ) + \text{i} \beta \left ( 2 \tilde {v}_{10} u_{20}-\tilde {u}_{10} v_{20} \right ) \biggr . \nonumber \\ & \quad \left . +\, \frac {\text{d} \tilde {u}_{10}}{\text{d} z} w_{20} + \tilde {w}_{10} \frac {\text{d} u_{20}}{\text{d} z} + w_{10} \frac {\text{d} U_1}{\text{d} z} \right ]\!, \end{align}

(4.24) \begin{align} \mathcal{G}_y & =- \textit{Re} \left [ \text{i} \alpha \left ( 2 \tilde {u}_{10} v_{20} - \tilde {v}_{10} u_{20} + v_{10} U_1 \right )+ \text{i} \beta \tilde {v}_{10} v_{20} + \frac {\text{d} \tilde {v}_{10}}{\text{d} z} w_{20} + \tilde {w}_{10} \frac {\text{d} v_{20}}{\text{d} z} \right ]\!, \end{align}

(4.25) \begin{align} \mathcal{G}_z & =- \textit{Re} \! \left [ \text{i} \alpha \left (2 \tilde {u}_{10} w_{20} - \tilde {w}_{10} u_{20} + w_{10} U_1 \right ) + \text{i} \beta \left ( 2 \tilde {v}_{10} w_{20}-\tilde {w}_{10} v_{20} \right ) + \frac {\text{d}}{\text{d} z} \left ( \tilde {w}_{10} w_{20} \right ) \! \right ]\!, \end{align}

where the symbol ‘ $\sim$ ’ in (4.5) and (4.23)–(4.25) denotes the complex conjugate. The consideration of terms above the orders of $c_i^{3/2}$ is not necessary for cubic Landau theory and, thus, is neglected in the present study (Khandelwal & Bera Reference Khandelwal and Bera2015; Khan, Bera & Khandelwal Reference Khan, Bera and Khandelwal2019).

The system of harmonic equations (4.5)–(4.22) is solved sequentially in the increasing powers of $c_i$ and $c_{ip}$ . Since the relations (3.3) and (4.4) indicate that the Fourier expansion coefficients in the fluid and porous region are related as $c^{1/2}_{i} B=c^{1/2}_{ip} B_{\!p}$ , therefore, to avoid the redundancy, the following discussion is given in terms of increasing powers of $c_i$ only.

At order $c_i^0$ , the $E^0$ harmonic equations are reduced to the basic state (2.7)–(2.8). Therefore, the terms $U_{0}$ and $U_{0p}$ denote the basic state velocities in fluid and porous regions, respectively. There are no contributions from $E^1$ and $E^2$ harmonics. At order $c_{i}^{1/2}$ , the $E^1$ harmonic contains the linear stability (3.4)–(3.11), and the contribution from $E^0$ and $E^2$ harmonics are zero. Therefore, the terms with ‘10’ subscript in the Fourier expansions (4.1) and (4.2) represent the solution obtained from the linear stability theory. At order $c_{i}^{1}$ , the equations corresponding to $E^0$ and $E^2$ harmonics contain the non-homogeneous system of linear equations to calculate the distorted mean flow and ‘20’ system, respectively. Since the non-homogeneous part on the right-hand side of these systems is known from the linear stability theory, the distorted mean function and the ‘20’ system can be easily computed. At order $c_{i}^{3/2}$ , the $E^1$ -harmonic (4.7)–(4.14) become non-homogeneous with the terms on the right-hand side containing the coefficients of ${\text{d}B}/{\text{d}\tau }$ , $B$ and $B | B |^2$ . These coefficients are known from the lower-order analysis. The corresponding homogeneous form of the above-mentioned system consists of the perturbed linear operator operating on the unknown functions with an ‘11’ subscript. For a non-trivial solution, the integrability condition requires that the non-homogeneous terms of $E^1$ -harmonic at order $c_i^{3/2}$ are orthogonal to the adjoint functions corresponding to the linearised disturbed system. The homogeneous adjoint equations of the linear system are given in the Appendix C. Therefore, multiplying the non-homogeneous parts in (4.7)–(4.14) with the corresponding adjoint variables $p^*, \, u^*, \, v^*, \, w^*, p_{\!p}^*, \, u_{\!p}^*, \, v_{\!p}^*, \, w_{\!p}^*$ , and using the normalisation condition (C17) yields the cubic Landau equation

(4.26) \begin{equation} \frac {\text{d} B}{\text{d} \tau } = \alpha B + a_1 B | B |^2, \end{equation}

where the first Landau coefficient $a_1$ is given by

(4.27) \begin{equation} a_1=\frac {1}{\textit{Re}} \int _{0}^{1} \big ( \mathcal{G}_x u^* + \mathcal{G}_y v^* + \mathcal{G}_z w^* \big ) \, \text{d}z . \end{equation}

Note that the Landau equation (4.26) gives a modification to the exponential growth/decay of the disturbance predicted by the linear stability analysis.

4.1.2. Physical significance of first Landau constant and amplitude function

The physical interpretation of (4.26) can be determined by the concept of equilibrium amplitude as follows.

Substituting $B=c_{i}^{1/2} A$ , $\tau =c_i t$ in the Landau equation (4.26), we get

(4.28) \begin{equation} \frac {\text{d} |A|^2}{\text{d} t}= 2 \alpha c_i |A|^2 + 2 a_{1r} |A|^4, \end{equation}

where $a_{1r}=\text{Real} ( a_1 )$ and $A$ is called physical amplitude of the disturbance wave. The stationary amplitude, $A_s$ of (4.28), is obtained by satisfying the stationary condition ${\text{d} | A |}/{\text{d} t}=0$ . This gives the possibility of three stationary amplitude solutions given by

(4.29) \begin{equation} A_s=0\, , \,\quad A_s=\pm \sqrt {-\frac {\alpha c_i}{a_{1r}}}. \end{equation}

Note that $A_s$ is always a positive quantity; therefore, the negative solution in (4.29) is insignificant. The nature of the other two equilibria becomes clear once the analytical solution to (4.28) is obtained, which is presented as follows (Stuart Reference Stuart1960; Khandelwal & Bera Reference Khandelwal and Bera2015):

(4.30) \begin{equation} |A|^2=\frac {A_0 c_i e^{2 \alpha c_i t}}{\left ( 1- {a_{1r}}/{\alpha } A_0 e^{2 \alpha c_i t}\right )}, \end{equation}

where $A_0$ is an arbitrary real constant.

Figure 9.

Variation of physical amplitude with time in the (a) supercritical regime and (b) subcritical regime.

For $c_i\gt 0$ , i.e. when $\textit{Re} \geqslant \textit{Re}_c$ , where $ \textit{Re}_c$ is the critical Reynolds numbers corresponding to least stable disturbance, the solution (4.30) will be meaningful for all $t$ only if $| A |^2\gt 0$ , which is possible when $a_{1r}\lt 0$ . In this scenario, $| A |^2 \rightarrow 0$ as $t \rightarrow -\infty$ , and tends an equilibrium value $-({\alpha c_i}/{a_{1r}})$ as $t \rightarrow \infty$ . The phase portrait diagram for the supercritical case is shown in figure 9(a), indicating that the solutions $A_s=0$ and $A_s=\sqrt {-({\alpha c_i}/{a_{1r}})}$ correspond to unstable and stable equilibria, respectively, when $c_i\lt 0$ and $a_{1r}\gt 0$ . Since any perturbation in the supercritical regime will grow or decay to an equilibrium value of $\sqrt {-({\alpha c_i}/{a_{1r}})}$ , the amplitude $A_s$ is called the equilibrium amplitude ( $A_e$ ).

For $c_i\lt 0$ , i.e. when $ \textit{Re} \leqslant \textit{Re}_c$ , the solution (4.30) will be meaningful for all $t$ when $a_{1r}\gt 0$ . In this case, $| A |^2 \rightarrow -({\alpha c_i}/{a_{1r}})$ as $t \rightarrow -\infty$ . At some finite time $t_0$ , if $|A|^2$ is slightly less than $-({\alpha c_i}/{a_{1r}})$ , the amplitude of the disturbance decays to zero. However, if $|A(t_0)|^2 \geqslant -({\alpha c_i}/{a_{1r}})$ , the disturbance grows in magnitude. The phase portrait diagram for the subcritical case is given in figure 9(b), which shows that the solutions $A_s=0$ and $A_s=\sqrt {-({\alpha c_i}/{a_{1r}})}$ correspond to stable and unstable equilibria, respectively. Since a minimum threshold amplitude of disturbance is required to cause the flow instability, the amplitude $\sqrt {-({\alpha c_i}{a_{1r}})}$ in the subcritical regime is called the threshold amplitude ( $A_t$ ).

4.4.1. Finite-amplitude behaviour along the instability boundary

In this subsection, we identify the nature of bifurcation by examining the real part of the Landau constant along the critical instability boundary with $L_x$ . Accordingly, the BP is the critical instability point (CIP) obtained from linear stability theory in figure 8.

Figure 10.

Impact of $L_x$ on first Landau constant $a_{1r}$ along the BP, when (a) $\hat{d}=0.08, \, 0.10, \, 0.12$ and (b) $\hat{d}=0.15, \, 0.20, \, 0.30$ .

Figure 10 illustrates the variation of $a_{1r}$ with the streamwise slip length ( $L_x$ ) along the instability curve at various depth ratios. Since the growth rate $c_i$ is zero along the critical curve, the nature of bifurcation is defined only on the basis of the sign of $a_{1r}$ . Accordingly, for a positive (negative) value of $a_{1r}$ , the bifurcation is termed as subcritical (supercritical). Note that, in the vicinity of $L_x=0$ , for $\hat{d}=0.08, \, 0.10$ (PM) and for $\hat{d}=0.12$ (OFM), the negative sign of $a_{1r}$ indicates the supercritical bifurcation. As the streamwise slip increases, the nature of bifurcation shifts from supercritical to subcritical in the regime of $L_x$ , where three-dimensional disturbances dominate the instability. Consequently, for $\hat{d}=0.08, \, 0.10$ and $0.12$ , $a_{1r}$ changes its sign from negative to positive at $L_x=0.035, \, 0.014$ and $0.004$ , respectively. A further increment in $L_x$ lowers the magnitude of $a_{1r}$ and thus leads back to the supercritical bifurcation. When the two-dimensional disturbances again become least stable, $a_{1r}$ increases monotonically with $\hat{d}$ . On increasing $\hat{d}$ beyond $0.12$ in figure 10(b), $a_{1r}$ is positive and increases monotonically with $L_x$ near zero. When three-dimensional disturbances become linearly least stable, $a_{1r}$ decreases steadily with $L_x$ , attains a local minimum within the supercritical bifurcation regime and then increases toward the positive values. As two-dimensional disturbances regain dominance, $a_{1r}$ remains positive throughout, leading to a subcritical bifurcation. Overall, we observe that the range of $L_x$ where the subcritical bifurcation exists is found to be increasing monotonically with $\hat{d}$ . Since the type of bifurcation is subcritical for the limiting channel case (Stewartson & Stuart Reference Stewartson and Stuart1971; Orszag & Patera Reference Orszag and Patera1980), such a behaviour is anticipated for high values of $\hat{d}$ . Note that the bifurcation characteristics along the instability curve for $\hat{d}$ equal to 0.08, 0.1, 0.12 are similar and, at the same time, the bifurcation characteristics for $\hat{d}$ equal to 0.15, 0.2, 0.3 are identical. Therefore, to understand the impact of the location of the interface on the nonlinear evolution, we have considered two values, 0.1 and 0.2 of $\hat{d}$ , for the rest of the analysis.

Figure 11.

Variation of nonlinear KE components with respect to $L_x$ along the BP, when (a) $\hat{d}=0.1$ and (b) $\hat{d}=0.2$ .

To uncover the driving mechanism behind the observed bifurcation characteristics, we employ the nonlinear kinetic energy (NKE) balance. Figures 11(a) and 11(b) illustrate the variation of different NKE components with $L_x$ along the instability curve for $\hat{d}=0.1$ and $0.2$ , respectively. In both figures, the sum of the components in the NKE spectra equals twice $a_{1r}$ , which supports the nonlinear stability results obtained in the present study under Landau theory. Furthermore, the magnitude of $K_{11}$ is almost negligible, indicating an insignificant change in the variation of KE via modification to the surface drag. At $L_x=0.014$ in figure 11(a), a drastic change in the NKE spectra is seen. Here, the magnitude of both $D_{11}$ and $P_{110}$ increases suddenly to large values and the positive sign of $D_{11}$ indicates a rapid decrease in the rate of dissipation of the disturbance. Apart from this, modification in viscous dissipation of disturbance kinetic energy due to a change in the shape of fundamental wave (i.e. $D_{11}$ ) surpasses the combined negative contributions from $P_{101}$ , $P_{110}$ and $E_{12}$ , resulting in the subcritical bifurcation. As $L_x$ increases beyond $0.014$ , the continuous drop in the value of $D_{11}$ causes a monotonic decrease in $a_{1r}$ and results in supercritical bifurcation for $L_x \geqslant 0.041$ . In figure 11(b), the positive contributions of $P_{110}$ and $D_{11}$ drive the subcritical bifurcation over the range $L_x \in [0, 0.008]$ . When $0.009 \leqslant L_x \leqslant 0.027$ , $2a_{1r}$ decreases continuously and attains its negative minima. Note that at this point, the contribution of $P_{110}$ in NKE is almost balanced by $D_{11}$ , and thus the nature of bifurcation is decided mainly by gradient production due to interaction between fundamental mode and the distorted mean flow strain (i.e. $P_{101}$ ). As $L_x$ increases further, $P_{110}$ becomes more favourable in shear production and takes the lion’s role in deciding the nature of bifurcation, thus causing the subcritical bifurcation. Additionally, for $L_x\geqslant 0.05$ , the bifurcation phenomena mainly depend upon $P_{110}$ , while the other components are stabilising in nature.

4.4.2. Bifurcation as a function of $\alpha$

It is quite possible that the instability, which is supercritical (subcritical) at some wavenumber, may become subcritical (supercritical) at another wavenumber. In such situations, any potentially unstable wave may grow and interact with other modes. So, investigation of the bifurcation phenomenon along the NSC is essential. To address this objective, we have considered three different values $0$ , $0.02$ and $0.05$ of $L_x$ at $\hat{d}=0.1$ . The NSCs in the $(\alpha , \textit{Re})$ -plane and the variation of $a_{1r}$ along these curves are presented in figures 12(a), 12(b) and 12(c) for $\{L_x, \beta _c\}$ equal to $\{0,0 \}$ , $\{0.02,1.35\}$ and $\{0.05,0.98\}$ , respectively. The solid (dashed) curve represents the variation of $\textit{Re}$ ( $a_{1r}$ ) along the neutral curve with dot representing the $\textit{Re}_c$ (corresponding $a_{1r}$ ). Under the no-slip condition, throughout the neutral curve, the negative value of $a_{1r}$ illustrates the supercritical bifurcation. Whereas, for $L_x=0.02$ , the subcritical bifurcation exists up to $\alpha \leqslant 1.94$ (that includes CIP) and beyond which the supercritical bifurcation persists. On further increment of $L_x$ to $0.05$ , supercritical bifurcation exists at CIP and shifts to subcritical for the relatively small values of wavenumber (see inset of figure 12 c). Therefore, the nature of bifurcation changes as one moves along the NSC.

Figure 12.

Neutral stability curves for ${\textit{Da}}=10^{-6}$ at $\hat{d}=0.1$ for the critical spanwise disturbance when (a) $L_x=0$ , (b) $L_x=0.02$ , (c) $L_x=0.05$ . The corresponding dashed curves represent the variation of $a_{1r}$ along the neutral curve.

4.4.3. Finite-amplitude behaviour as a function of $\textit{Re}$

The cubic Landau theory is valid only in the vicinity of BP and traces the evolution of finite-amplitude perturbations. However, according to Rogers, Moulic & Yao (Reference Rogers, Moulic and Yao1993), the analysis can also be extended beyond the BP in the limit as $c_i$ approaches zero. It has been reported that in the linearly stable regime, the results obtained by the cubic landau theory, quintic Landau theory and DNS are qualitatively as well as quantitatively the same for ${O}(c_i)\lt 10^{-3}$ , and only qualitatively for ${O}(c_i)\lt 10^{-2}$ . However, in the linearly unstable region, the prediction from cubic Landau theory, even for $c_i \leqslant 4 \times 10^{-2}$ , is well matched with the DNS results.

Being motivated by these observed validity limits, we now impose an upper bound on $c_i$ in the present analysis for examining the bifurcation characteristics away from the CIP. We have seen from our numerical experiments that in the linearly stable regime, each $a_{1r}$ curve experiences a sudden jump at some $\textit{Re}$ , where the magnitude of $a_{1r}$ becomes almost zero and does not change as one moves further away from the CIP. If we define the corresponding value of $c_i$ as $c_i^*$ , then it has been checked that the maximum absolute value of $c_i^*$ for the parameters considered in this study is of the order of $10^{-2}$ and the $a_{1r}$ curves are traced as a function of $\textit{Re}$ accordingly. Note that a complete picture of nonlinear saturation can be confirmed by the use of direct numerical simulation or higher order nonlinear analysis, which is beyond the scope of the present work. Hence, the following analysis will provide a qualitative characterisation of the nonlinear interaction between different harmonic modes for larger (lower) values of $\textit{Re}$ beyond (before) the bifurcation point.

The nonlinear characteristics are analysed as a function of $\textit{Re}$ in terms of a parameter, $\delta _{\textit{Re}}$ , defined by $\delta _{\textit{Re}}= ({\textit{Re}}/{\textit{Re}_c})-1$ . Accordingly, $\delta _{\textit{Re}}\lt 0$ signifies the linearly stable regime, while $\delta _{\textit{Re}}\gt 0$ indicates the linearly unstable regime. The results are shown for $0.1$ and $0.2$ of $\hat{d}$ . Figure 13 shows the variation of $a_{1r}$ and stationary amplitude ( $A_s$ ) as a function of $\delta _{\textit{Re}}$ for different values $(0, 0.01, 0.03, 0.05, 0.08)$ of $L_x$ . Figures 13(a) and 13(b) correspond to $\hat{d} = 0.1$ , while figures 13(c) and 13(d) correspond to $\hat{d} = 0.2$ . We see for both values of $\hat{d}$ in the linearly stable zone, each $a_{1r}$ curve experiences a point of discontinuity and the corresponding value of $c_i^*$ is given in table 2. It can be seen from figure 13(a) that for the no slip as well as slip with $L_x = 0.01$ cases, the sign of $a_{1r}$ is negative in the entire range of $\delta _{\textit{Re}}$ , indicating the existence of supercritical bifurcation in the linearly unstable regime which introduces stable supercritical branch of equilibrium amplitude (see figure 13 b). This also indicates the existence of a new laminar flow for $\textit{Re} \gt \textit{Re}_c$ . For all other values of $L_x$ , the sign of $a_{1r}$ changes from positive to negative on moving from the left-hand side to the right-hand side of the bifurcation point. It shows that both subcritical and supercritical bifurcations are possible. Therefore, branches of both stable supercritical equilibrium and unstable subcritical threshold amplitudes can be seen in figure 13(b). For $L_x$ equal to 0.05 and 0.08, an unstable subcritical (stable supercritical) branch appears away from (at) the bifurcation point. Interestingly, at $L_x=0.03$ , the sign of $a_{1r}$ is positive in the range $-0.386 \leqslant \delta _{\textit{Re}} \leqslant 0.243$ ; consequently, the unstable branch of threshold amplitude appears in that range starting from the CIP and a stable branch of equilibrium amplitude appears in the unstable regime beyond $\delta _{\textit{Re}} = 0.243$ . As the location of the interface is shifted downward by changing the value 0.1 of $\hat{d}$ to 0.2, except for $L_x=0.03$ , the sign of $a_{1r}$ remains positive in the linearly stable regime and the threshold amplitude exists starting from the CIP (see figures 13 c and 13 d). For $L_x=0.03$ , the same starts from a point in the neighbourhood of CIP.

Table 2.

Critical point, $c_i^*$ and $\textit{Re}^*_{_{lc}}$ for different values of $\hat{d}$ and $L_x$ . The bold value represents when the BP for two-dimensional disturbances give a smaller $\textit{Re}^*_{_{lc}}$ as compared with the same at the CIP.

Figure 13.

Variation of $a_{1r}$ and $A_s$ with $\delta _\textit{{Re}}$ for different values of $L_x$ , when (a,b) $\hat{d}=0.1$ and (c,d) $\hat{d}=0.20$ .

Figure 14.

Variation of nonlinear KE components with respect to $\delta _{\textit{Re}}$ , when $\hat{d}=0.1$ for (a) $L_x=0$ , (b) $L_x=0.03$ , (c) $L_x=0.08$ ; and when $\hat{d}=0.2$ for (d) $L_x=0$ , (e) $L_x=0.03$ , ( f) $L_x=0.08$ .

Thus, we see that the impact of slip on the finite-amplitude solution in the neighbourhood of BP is highly sensitive to the location of the interface. A nonlinear kinetic energy spectrum is used to understand the physics underlying this sensitivity characteristic. For this, we have chosen three values 0, 0.03, 0.08 of $L_x$ for both values 0.1 and 0.2 of $\hat{d}$ in figures 14(a–c) and 14(d–f), respectively. As seen in § 3.3.3, the least stable mode at $L_x$ equal to 0, 0.03, 0.08 are 2-D, 3-D and 2-D, respectively. In figure 14(a), when no slip condition is used, except for $D_{11}$ and $K_{11}$ , all other terms in the NKE spectrum are negative and $P_{101}$ is most dominating. Thus, even though $D_{11}$ and $K_{11}$ are not in favour of dissipating the disturbance, the energy released by the fundamental disturbance to distort the mean flow through $P_{101}$ is so prominent that the overall sum is negative and shows supercritical bifurcation in the regime of $\delta _{\textit{Re}}$ . For $L_x$ equal to $0.03$ , $P_{110}$ ( $D_{11}$ ) becomes the most stabilising (destabilising) factor in the entire regime and its magnitude increases as one moves to the right-hand side of BP (see figure 14 b). The nature of subcritical or supercritical bifurcation depends on whether a decrease in the rate of dissipation of disturbance dominates over the strong unwillingness in shear production of $P_{110}$ or vice versa. In the case of $L_x$ equal to 0.08 in figure 14(c), the nature of $P_{110}$ ( $D_{11}$ ) is reversed. Consequently, subcritical phenomena are observed away from the BP, whereas supercritical phenomena are seen starting from the BP. In contrast, at $\hat{d} = 0.2$ , $P_{110}$ is positive and increases monotonically in entire permissible range of $\delta _{\textit{Re}}$ for $L_x$ equal to $0$ and $0.08$ . Furthermore, the destabilising nature of $P_{110}$ also increases as the interface shifts downward. For $L_x=0.03$ in figure 14(e), $P_{110}$ and $D_{11}$ are decreasing in nature. Consequently, the nature of bifurcation as a function of $\delta _{\textit{Re}}$ resembles that of the case of $L_x=0.08$ when $\hat{d}=0.1$ . Overall, we observe that the characteristics of $D_{11}$ and $P_{110}$ decide the nature of the bifurcation. It has also been checked through our numerical experiments that as we further increase the value of $\hat{d}$ , the contribution of $P_{110}$ becomes increasingly more prominent and the nature of bifurcation remains subcritical, irrespective of the value of slip length.

Before discussing the role of finite-amplitude evolution on physical quantities, such as skin-friction, in the next subsection, an important note is made regarding the existence of a lower critical point. According to Landau’s suggestion, a lower critical point may exist in the subcritical regime. Note that the lower critical point is a threshold value of $\delta _{\textit{Re}}$ at which the unstable branch turns into a stable branch, which is the saddle node point, and below it, the bifurcated solution does not exist. To find such a point, one has to rely either on the DNS or higher-order quintic-Landau theory. The way the present analysis is made, the quintic Landau equation will have three stationary solutions, out of which two are stable and one is unstable. The point where the non-zero stable and unstable branches meet is called a saddle node point, and the corresponding $\delta _{\textit{Re}}$ gives the lower critical $\textit{Re}$ . The analysis based on the fifth-order Landau equation is beyond the scope of the present study. However, by the way $a_{1r}$ changes as a function of $\textit{Re}$ in the present study, it is clear that the lower critical value of $\textit{Re}$ cannot be smaller than $\textit{Re}$ calculated at $c_i^*$ , say $\textit{Re}_{lc}^*$ . Thus, $\textit{Re}_{lc}^*$ is defined as the probable lower critical corresponding to $\textit{Re}_c$ . The critical points, $c_i^*$ , and the corresponding $\textit{Re}_{lc}^*$ for different set values of $(\hat{d}, L_x)$ are illustrated in table 2. Interestingly, for $\hat{d}=0.2$ and $L_x = 0.01$ , it has been found that although the critical $\textit{Re}$ under three-dimensional disturbance is smaller than the same under two-dimensional disturbance, $\textit{Re}_{lc}^*$ obtained under two-dimensional disturbance is lower than the same obtained under three-dimensional disturbance (see the bold entry in table 2).

Overall, the weakly nonlinear stability analysis reveals a significant relationship between the critical mode of instability and the nature of bifurcation. The PM and OFM are always characterised by the supercritical bifurcation. However, for the EFM, two scenarios are possible: (i) dual existence of subcritical and supercritical bifurcation; or (ii) only subcritical bifurcation, depending upon the set of parameters involved. The presence of a lower critical point in the subcritical regime indicates an early onset of instability in the flow.

4.4.4. Variation of skin-friction coefficient

The transition from laminar to turbulent is always accompanied by a change in the skin-friction coefficient. Various experimental studies available in the literature have been performed to calculate the critical value of controlling parameters, such as Reynolds number, at which a change in friction drag occurs. In this regard, we now analyse the variation of the skin-friction coefficient for the basic and distorted state under the impact of finite-amplitude disturbances. The skin-friction coefficients at the upper wall and the interface are given by (Sharma et al. Reference Sharma, Khandelwal and Bera2018)

(4.41) \begin{equation} C_{w}=-\frac {2}{\textit{Re}} \left .\frac {\text{d}U}{\text{d}z}\right |_{z=1}, \quad C_{i}=\frac {2}{\textit{Re}} \left .\frac {\text{d}U}{\text{d}z}\right |_{z=0}, \end{equation}

respectively. The velocity function $U$ for the distorted state is determined at the stationary amplitude, i.e. $U_0+ A_{s}^2 U_1$ . Figure 15 illustrates the impact of streamwise slip on the basic and distorted state for the skin-friction factor at the wall and the interface. In figure 15(a), the role of $L_x$ on friction coefficient for the basic state at the interface and wall (i.e. $C_{ib}$ and $C_{wb}$ , respectively) is analysed for two different values, 0.1 and 0.2, of $\hat{d}$ at constant Reynolds number of $5000$ . Irrespective of the value of the slip coefficient, an increase in $\hat{d}$ causes an increment in the value of the skin-friction coefficient at both the wall and the interface, while an increase in $L_x$ causes a simultaneous decrement in the value of $C_{wb}$ and $C_{ib}$ . Since the slip is imposed on the upper wall, the change in friction coefficient with $L_x$ at the upper wall is steeper than the same at the interface. Consequently, the depth ratio increases the skin-friction coefficient, whereas the slip decreases the same at both the wall and the interface. To see the impact of finite amplitude disturbances on the skin-friction coefficient at the top plate and the interface, and its dependency on slip length, figures 15(b) and 15(c), respectively, have been plotted for different values of $L_x$ . It is to be noted that the friction factor is not defined in the vicinity of $\delta _{\textit{Re}}$ where $a_{1r}$ changes its sign. At this point, the finite-amplitude solution does not exist. Irrespective of the slip coefficient, figure 15(b) shows that the skin-friction for the distorted state is higher than that for the base state. However, the skin-friction at the wall (i.e. $C_{wd}$ ) may be higher or lower than the skin-friction calculated under the base state (i.e. $C_{wb}$ ). In general, it can be seen that when the linearly least stable disturbances are two-dimensional (three-dimensional), $C_{wd}$ is always less (greater) than $C_{wb}$ . Quantitatively, we have observed that the increase (decrease) in the skin friction at the interface (upper wall) due to distortion of base flow is around 15.1 % (4.5 %) for $L_x=0.08$ at $\delta _{\textit{Re}}=0.4$ .

Figure 15.

Variation of skin-friction coefficient with $L_x$ on the (a) basic state at the wall and interface for $\hat{d}=0.1$ and $0.2$ when $\textit{Re}=5000$ , (b) basic and distorted state at the fluid–porous interface for $\hat{d}=0.1$ as a function of $\delta _{\textit{Re}}$ , and (c) basic and distorted state at the wall as a function of $\delta _{\textit{Re}}$ . In panels (b) and (c), the dashed lines are skin-friction coefficients for the distorted flow.

Therefore, the transition is observed by an increase in the skin-friction coefficient at both the upper wall and interface when three-dimensional disturbances dominate the linear stability, whereas the same is found by an increase at the interface only when two-dimensional disturbances are linearly least stable. Our numerical experiments reveal that similar observations are found for other values of $\hat{d}$ .