Turbulent transition in a channel with superhydrophobic walls: anisotropic slip and shear misalignment effects

Abstract

Superhydrophobic surfaces dramatically reduce skin friction of overlying liquid flows. These surfaces are complex and numerical simulations usually rely on models to reduce this complexity. One of the simplest consists of finding an equivalent boundary condition through a homogenisation procedure, which in the case of channel flow over oriented riblets, leads to the presence of a small spanwise component in the homogenised base flow velocity. This work aims at investigating the influence of such a three-dimensionality of the base flow on stability and transition in a channel with walls covered by oriented riblets. Linear stability for this base flow is investigated: a new instability region, linked to cross-flow effects, is observed. Tollmien–Schlichting waves are also retrieved but the most unstable are three-dimensional. Transient growth is also affected as oblique streaks with non-zero streamwise wavenumber become the most amplified perturbations. When transition is induced by Tollmien–Schlichting waves, after an initial exponential growth regime, streaky structures with large spanwise wavenumber rapidly arise. Modal mechanisms appear to play a leading role in the development of these structures and a secondary stability analysis is performed to retrieve successfully some of their characteristics. The second scenario, initiated with cross-flow vortices, displays a strong influence of nonlinearities. The flow develops into large quasi-spanwise-invariant structures before breaking down to turbulence. Secondary stability on the saturated cross-flow vortices sheds light on this stage of transition. In both cases, cross-flow effects dominate the flow dynamics, suggesting the need to consider the anisotropicity of the wall condition when modelling superhydrophobic surfaces.

1. Introduction

Drag is the resistance to motion experienced by a fluid flowing on a surface, generated by the difference in velocity between the solid object and the fluid. The reduction of drag represents a key factor for a whole range of physical and engineering problems involving the relative motion of a fluid on a solid surface, for instance the transport of drinkable water in pipes, of the blood in human vessels, and the motion of aircraft in the air and of ships in the sea. For all of these applications, one of the main sources of drag is the skin friction between the molecules of the fluid and the solid surface over which it flows, whose magnitude depends on the properties of both the surface and the fluid flowing on it. In the vicinity of a smooth hydrophilic surface, the flow must decelerate until reaching zero velocity at the wall, self-inducing a strong resistance to motion due to skin friction. This resistance to motion can be reduced by decreasing the gradient of velocity between the surface and the flow itself, which can be achieved by using particular surfaces or coatings, able to assure shear-free or slip wall conditions. One example of this kind of surfaces is given by superhydrophobic (SH) surfaces (Rothstein 2010). The superhydrophobicity of a surface is due to its nanostructure, which is composed by a hierarchical structure of microroughnesses that trap the air underneath them, reducing the surface of contact with the water droplets and the wetting of the surface. Such a hierarchical, rough nanostructure can be found on many biological surfaces such as those of lotus leaves, butterfly wings, duck feathers and water striders legs, and have inspired the engineering of biomimetic non-wettable materials for applications that range from self-cleaning to anti-icing.

The experimental works of Tretheway & Meinhart (2002) and Choi, Westin & Breuer (2003) were among the first to observe the drag-reducing effect of a SH surface. These studies, focusing on laminar flows, demonstrated the existence of large slip velocities in the vicinity of the gas pockets. In view of these promising results, several other experimental studies were performed to better understand the underlying mechanisms. Since drag-reducing properties appeared to depend strongly on the geometry of the SH surface, one of the first objectives was to quantify (Truesdell et al. 2006; Tsai et al. 2009) or even maximise (Lee & Kim 2009) the slip velocities and associated drag reduction for a given configuration. For example, using different types of micro-posts and micro-ridges, Ou, Perot & Rothstein (2004) have been able to show that an increasing shear-free area induced by carefully engineered micro-roughnesses induces an increased slip length and a consequent drag reduction for laminar channel flows. Along the air–water interfaces, the authors measured slip velocities larger than the $60\,\%$ of the averaged velocity, corresponding to a drag reduction larger than $40\,\%$. More experiments have been performed using SH surfaces characterised by nano-sized structures such as ridges or needles (Choi & Kim 2006), trying to engineer the nano-structures able to maximise the attainable drag reduction (Lee & Kim 2009). In particular, Steinberger et al. (2007) have found that microposts are less effective in reducing drag than ridges, since the flow has to decelerate and accelerate between different posts resulting in a lower slip length. Direct numerical simulations (DNS) have followed these experiments; several levels of complexity were progressively reached. Initially, DNS were performed by imposing a slip condition using an arbitrary slip length (Min & Kim 2004). Later, in order to model the interaction between air pockets and water, the air–water interface was assumed flat and the viscosity of the air trapped in the micro-ridges was neglected (Ybert et al. 2007). The resulting boundary condition consisted in an alternation between shear-free and no-slip patches on a flat surface (Martell, Perot & Rothstein 2009; Martell, Rothstein & Perot 2010). More recent simulations take interface deformation into account, either modelled through the use of a linearised Young–Laplace equation (Seo, García-Mayoral & Mani 2018), either directly by simulating the trapped lubricant altogether (Bernardini et al. 2021; Sundin, Zaleski & Bagheri 2021).

Comparisons of these different approaches have shown that, as long as the texture size remains small, using spatially homogeneous partial slip boundary conditions is a reliable approach in SH surfaces modelling: it is capable of correctly predict properties for both laminar (Choi et al. 2003; Ou et al. 2004) and turbulent regimes (Zhang et al. 2015; Zhang, Yao & Hao 2016). Stability thresholds and instability mechanisms are also accurately retrieved (Yu, Teo & Khoo 2016). Still, the accurate description of SH surfaces rests on solving two fundamental difficulties: finding an adequate slip length modelling the effect of the SH wall and ensuring that this slip length describes a wetting-stable configuration within the micro-roughnesses (namely, one able to retain the trapped lubricant even in turbulent conditions). The stability of an air–water interface in the context of SH coatings is complex and based on several phenomena, as studied parametrically by Seo et al. (2018). Wetting transition, that causes plastron depletion and ultimately leads to augmented drag, may be caused by capillarity or stagnation pressure effects. Recently, Seo et al. (2018) have shown that for limited values of the Weber number and of the roughness size, wetting stable conditions are maintained, while choosing larger values of these parameters may lead to destabilisation of the liquid–gas interface. Considering the value of the slip length required for accurately modelling a SH surface, once a type of roughness and a set of parameters are chosen ensuring wetting-stable conditions, one can use a homogenisation approach. The first application of such an approach was proposed in the work of Sarkar & Prosperetti (1996) for a rough surface. Later, a very general framework, valid for any type of textured surface, was proposed by Bazant & Vinogradova (2008) in which the existence of a linear relation between the velocity and the shear stress components at the virtual interface is supposed, thus generalising the Navier slip boundary condition (Ybert et al. 2007). A method to explicitly compute the different slip lengths from the roughness geometry can be found in the article of Kamrin, Bazant & Stone (2010). Further progress on small-periodicity roughnesses was made in the work of Luchini (2013). Ultimately, higher-order relations were obtained in the works of Bottaro (2019), Bottaro & Naqvi (2020) and Sudhakar et al. (2021).

Park, Park & Kim (2013) have found that the potential drag reduction of SH coatings is much larger in turbulent flows than in laminar flows, and this effect might be due to the damping of wall turbulence induced by the presence of a slip length. However, the effect of those kind of surfaces on turbulent transition is a point which has been still not widely investigated. In general, transition from laminar to turbulent flow induces a strong increase in skin friction, along with a strong increase of the drag. Thus, in the transitional regime, the competition between drag decrease due to the surface micro-structure, and drag increase due to transition to turbulence might provide surprising results. Up to now, very few studies have been performed on transition to turbulence of a laminar flow over SH surfaces, focusing at first on the very first phase of transition, namely, linear instability (Lauga & Cossu 2005; Min & Kim 2005). For a channel flow, it has been proved by a local instability analysis that, when imposing a simple slip condition, the onset of two-dimensional (2-D) Tollmien–Schlichting (TS) waves is considerably postponed, allowing the flow to stay laminar up to larger Reynolds numbers, and further decreasing the drag. However, shear flows very often experience subcritical transition to turbulence, due to the transient growth of non-modal disturbances, bypassing the asymptotic growth of TS waves (Schmid & Henningson 2001; Farano et al. 2016). In particular, it has been shown in Min & Kim (2005) that slip boundary conditions have a strong influence on the linear growth of modal disturbances, but a very weak influence on the maximum transient energy growth of perturbations at subcritical Reynolds numbers, concluding that slip boundary conditions are not likely to have a significant effect on the transition to turbulence in channel flows. By means of global stability analyses, Tomlinson & Papageorgiou (2022) have shown that in the presence of SH grooves, additional instabilities may arise, with critical Reynolds numbers small enough to be achievable in applications.

The validity of linear stability results on laminar–turbulent transition itself, which is an intrinsically nonlinear phenomenon, has been recently assessed by Picella, Robinet & Cherubini (2019, 2020), who have confirmed that SH surfaces strongly influence transition induced by wall-close disturbances, such as TS waves, even at subcritical Reynolds number, but have a weak effect on the subcritical growth of coherent structures lying farther from the wall, such as streaks and streamwise vortices. Cherubini, Picella & Robinet (2021) reported a strong effect of boundary slip on the transient growth of nonlinear optimal perturbations: in particular, while the maximal energy growth is considerably decreased, the shape of the optimal perturbation barely changes, indicating the robustness of optimal perturbations with respect to wall slip.

Most literature studies dealing with the effect of SH surfaces on turbulent transition consider the case of statistically isotropic micro-roughnesses, therefore allowing for the use of a single slip length, $L_s$, acting on both wall-parallel directions. In a turbulent setting, exceptions include the works of Min & Kim (2004), Aghdam & Ricco (2016) and, more importantly, Busse & Sandham (2012) who considered anisotropic slip boundary conditions. Ultimately, the influence of slip lengths, either streamwise or spanwise, on the overlying turbulence remains limited (Ibrahim et al. 2021). From a transitional point of view, the influence of the homogenised boundary condition was initially investigated by Min & Kim (2005) who found that spanwise slip induces earlier transition while streamwise slip considerably delays it. More recently, Chai & Song (2019), Xiong & Tao (2020) and Zhai, Chen & Song (2023) demonstrated the link between the presence of spanwise slip and the instability of a three-dimensional (3-D) TS wave. However, these studies considered one value of the slip length for both the streamwise and spanwise directions, which was then nullified on one of the two wall-parallel directions, while not considering the case of different (non-zero) values of the slip. Conversely, in the case of oriented (non-isotropic) micro-posts, the slip tensor is no longer aligned with the main flow direction. Homogenised boundary conditions are then anisotropic and induce a shear misalignment in the vicinity of the wall. In addition to linear stability, the effect on turbulent transition of such tensorial slip conditions and the subsequent shear misalignment, which are representative of non-isotropic SH micro-roughnesses, remains entirely to be unveiled.

In order to investigate this issue, in the present paper we consider the case of oriented SH riblet-like micro-roughnesses such as those considered in Pralits, Alinovi & Bottaro (2017), and investigate their effect on modal transition. We show that the non-isotropicity of the surface gives origin to a new modal instability, similar to cross-flow instabilities recovered on the flow over swept wings, and will considerably alter non-modal stability as well. The transition to turbulence originated by the former instability, as well as by the slip-modified TS waves, are studied in detail using DNS and secondary stability analyses.

The paper is structured as follows. In § 2 we present the governing equations and the methods used to implement the SH surface, namely homogenised boundary conditions. In § 3 we show, using stability and transient growth analysis, how the dynamics of infinitesimal perturbations is influenced by the anisotropy introduced by the boundary conditions. DNS are performed to better understand the mechanisms at play: § 4 reports the numerical configuration of the case. In § 5, the DNS is initialised with a 3-D TS wave while a cross-flow like mode is used in § 6. In both sections, secondary stability analyses are carried out, shedding some light on the influence of cross-flow components in the transition processes. A final discussion and conclusions are given in § 7.

2. Governing equations

The flow of an incompressible Newtonian fluid in a channel of height $2h$ with SH walls is considered. The reference frame is chosen as $({x},{y},{z})$, with $x$ being aligned with the pressure gradient, $y$ the wall normal direction and $z$ the direction orthogonal to the pressure gradient. In the following, for brevity, $x$ and $z$ are respectively denoted as the streamwise and spanwise directions.

For non-dimensionalising the problem, we should choose reference quantities, such as the half channel height $h$ for the length, a reference velocity $U_r$ and the kinematic viscosity of the fluid $\nu$, so that the Reynolds number is defined as $Re = U_rh/\nu$. Concerning the reference velocity, following previous works (Pralits et al. 2017; Chai & Song 2019; Xiong & Tao 2020; Zhai et al. 2023) we decided to set $U_r = 3U_a/2$, with $U_a$ being the average of the base flow velocity ($U_a = 1/2h \int U \,{{\rm d}y}$), so that the Reynolds number remains fixed at constant flow rate. Note that $U_r$ corresponds to the maximum velocity in the flow field in the laminar smooth-wall case, whereas in the SH cases it does not represent any intrinsic property of the flow, so that it is chosen as reference velocity mostly for comparison issues. After non-dimensionalisation, we obtain the following Navier–Stokes (NS) equations governing the behaviour of the flow:

(2.1)$$\begin{gather} \frac {\partial{\boldsymbol{U}}} { \partial{t}} ={-} (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U} -\boldsymbol{\nabla} p +\frac{1}{Re}{\nabla}^2 \boldsymbol{U} , \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}$$

where $\boldsymbol {U}=(U,V,W)^{\rm T}$ is the velocity vector and $p$ is the pressure. The flow is periodic in the streamwise and spanwise directions, and different domain sizes will be considered, as detailed in § 4. The walls of the channel are covered with SH riblet-like roughnesses oriented with an angle $\theta$, defined with respect to the $x$ direction (see figure 1(a)).

Figure 1. (a) Sketch of the SH surface covered with riblets oriented at an angle $\theta$ with respect to the streamwise direction. (b) Streamwise (orange) and spanwise (blue) velocity components of the base flow for $\lambda _{\parallel } = 0.03$ and $\theta =45^\circ$. (c) Angle between the bulk velocity and the pressure gradient.

For $\theta = 0^\circ$, the grooves are longitudinal whereas $\theta =90^\circ$ corresponds to transverse riblets. The effect of longitudinal SH riblet-like roughnesses can be modelled using equivalent streamwise and spanwise slip lengths $\lambda _{\parallel }$ and $\lambda _{\perp }$ (Gogte et al. 2005; Belyaev & Vinogradova 2010), which lead to the homogenised boundary conditions $U=\lambda _{\parallel }\partial _y U$, $V=0$ and $W=\lambda _{\perp }\partial _y W$. When the grooves are aligned with the mean pressure gradient, the streamwise slip length is twice the spanwise one, i.e. $\lambda _{\parallel } = 2 \lambda _{\perp }$ (Philip 1972; Asmolov & Vinogradova 2012). On the other hand, when the riblets are not aligned with the pressure gradient, the homogenised boundary conditions for the NS equations have the following more general form (Bazant & Vinogradova 2008),

(2.3)\begin{equation} \begin{bmatrix} U \\ W \end{bmatrix} =\mp \boldsymbol{L} \partial_y \begin{bmatrix} U \\ W \end{bmatrix}, \quad V = 0, \quad \text{at} \ y=\pm 1 \end{equation}

where the mobility tensor $\boldsymbol {L}$ depends on $\lambda _{\parallel }$ and on the rotation matrix $\boldsymbol {R}(\theta )$, allowing the rotation of the surface of an angle $\theta$:

(2.4)\begin{equation} \boldsymbol{L} = \boldsymbol{R} \begin{pmatrix} \lambda_{||} & 0 \\ 0 & \lambda_{{\perp}} \end{pmatrix} \boldsymbol{R}^T= \frac{\lambda_{{\parallel}}}{2} \begin{pmatrix} 1 + \cos^2 \theta & \cos{\theta}\sin\theta \\ \cos{\theta}\sin\theta & 1 + \sin^2 \theta. \end{pmatrix} \end{equation}

Following Pralits et al. (2017), we set $\lambda _{\parallel }=0.03$ and $\theta =45^\circ$ yielding $L_{11} = L_{22} = 0.0225$ and $L_{12} = L_{21} = 0.0075$. The choice of the slip length is critical and needs to be discussed. For linear stability analyses, this point has been partially assessed in the study of Yu et al. (2016) in which two modelling approaches of the SH surface have been used. In the first case, the SH surface has been replaced by an isotropic slip boundary condition while in the second configuration, the boundary conditions consist of an alternation of slip and no-slip patches. Note that the interface dynamics is neglected in this second case. This may be justified as DNS performed by Picella et al. (2019) showed that, for a supercritical transition, it played a limited role. Introducing the texture periodicity $L$ and the shear-free fraction $\delta$ of the surface, it was found that, for moderate periodicities ($L\geq 0.1$) and almost all $\delta$, an interfacial mode not captured with the homogenised approach arises and exhibits a much stronger growth rate than the typical TS waves. For periodicity an order of magnitude below ($L \approx 0.02$), such a mode does not develop no matter the value of the shear-free fraction and a good agreement between the two modelling approaches is retrieved. From there, the equivalent slip length of the SH surface can be expressed for small $L$ as a function of $L$ and $\delta$ only (Lauga & Stone 2003):

(2.5)\begin{equation} \lambda_{{\parallel}} = \frac{2L}{\rm \pi}\ln\left(\sec\left(\frac{\delta {\rm \pi}}{2}\right)\right). \end{equation}

Based on the previous discussion, setting $L = 0.025$ and $\delta = 0.90$ in the previous equation gives an equivalent slip length $\lambda _{\parallel }=0.03$. The configuration chosen in this article thus corresponds to a reasonable case in which the texture periodicity remains small enough to ensure the validity of the homogenised approach for a linear stability analysis.

Assuming a constant pressure gradient in the streamwise direction, the laminar steady base flow takes the form $\boldsymbol {U}_0 = (U_0(y),V_0,W_0)$ and is governed by the following equations:

(2.6a–c)\begin{equation} \frac{1}{Re}\frac{{\rm d}^2U_0}{{{\rm d}y}^2} ={-}\frac{{\rm d}P}{{\rm d}\,x}, \quad V_0 = 0, \quad \frac{{\rm d}^2 W_0}{{{\rm d}y}^2} = 0. \end{equation}

Solving the system together with the homogenised boundary conditions (2.3), the base flow velocity components read

(2.7a,b)\begin{equation} U_0 = \frac{ 2L_{11} + 1 - y^2}{1 + 3L_{11}}, \quad W_0 = \frac{2L_{12}}{1 + 3L_{11}}. \end{equation}

Figure 1(b)–(c) depicts the velocity components of the base flow and, more interestingly, the angle between the bulk velocity and the pressure gradient for $\lambda _{\parallel }=0.03$ and $\theta =45^\circ$. The effect of the spanwise component of velocity appears quite clearly: in the vicinity of the walls, the flow has a non-negligible angle with respect to the pressure gradient whereas, in the centre of the channel, the deviation caused by the presence of $W_0$ is less than $1^\circ$. The streamwise component of the base flow is similar to that found for a Poiseuille flow with a partially slippery wall with $\lambda _{\parallel } = \lambda _{\perp }$ (Philip 1972; Picella et al. 2019). Note that $U_0$ has the same mass flux as a plane Poiseuille flow with no slip, since for definition the Reynolds number remains fixed when the slip length is changed. The presence of a spanwise base flow component, which remains constant across the wall-normal direction may seem counterintuitive. However, being one order of magnitude smaller than the streamwise component, the spanwise one will mostly affect the near-wall regions, while having negligible effects in the centre of the channel.

3. Linear stability analysis

The instantaneous flow field is now decomposed as the sum of the previously described steady base flow and an unsteady disturbance having small amplitude. Introducing the state vector $\boldsymbol {Q}(\boldsymbol {x},t) = [U,V,W,P]^{\rm T}$, we have $\boldsymbol {Q}(\boldsymbol {x},t) = \boldsymbol {Q}_0(\boldsymbol {x}) + \epsilon \boldsymbol {q}(\boldsymbol {x},t)$, with $\epsilon \ll 1$. NS equations are then linearised with respect to the base flow $\boldsymbol {Q}_0=[\boldsymbol {U}_0,P_0]^{\rm T}$, yielding the following system of equations:

(3.1)$$\begin{gather} \frac {\partial{\boldsymbol{u}}} { \partial{t}} + (\boldsymbol{U}_0\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} + ( \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{U}_0 ={-}\boldsymbol{\nabla}p + \frac{1}{Re}\nabla^2 \boldsymbol{u}, \end{gather}$$

(3.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0 \end{gather}$$

together with the boundary conditions for the velocity perturbations,

(3.3)\begin{equation} \begin{bmatrix} u \\ w \end{bmatrix} =\mp \boldsymbol{L} \frac {\partial{}} { \partial{y}} \begin{bmatrix} u \\ w \end{bmatrix}, \quad v = 0 \quad \text{at} \ y={\pm} 1. \end{equation}

3.1. Modal stability analysis

The system being periodic in the streamwise and spanwise directions, the perturbation state vector $\boldsymbol {q}$ can be expanded in normal modes such that

(3.4)\begin{equation} \boldsymbol{q}(\boldsymbol{x},t) = \boldsymbol{\hat{q}}(y)\exp{[{\rm i}(\alpha x + \beta z - \omega t)]} + \text{c.c.} \end{equation}

with c.c. denoting the complex conjugate, $\alpha$, $\beta$ being the streamwise and spanwise wavenumbers, respectively, and $\omega$ an angular frequency. Temporal stability is investigated, implying $\alpha$ and $\beta$ are real whereas $\omega = \omega _r + {\rm i}\omega _i$ is a complex number whose imaginary part determines the asymptotic stability of the base flow $\boldsymbol {U}_0$ for a given mode $(\alpha, \beta )$. Thus, substituting (3.4) into the linearised NS equations (3.1)–(3.2) results in a generalised eigenvalue problem which can be solved by means of a spectral collocation code (Trefethen 2000; Schmid & Brandt 2014). The wall-normal direction is discretised with a Chebyshev grid. A grid constituted by $n=200$ points is used, ensuring a good convergence of the eigenmodes. The code is validated in the case of a smooth channel flow and compared with results from the book of Schmid & Henningson (2001).

In the present configuration, Squire's theorem does not hold. The demonstration for the present case can be seen in the article of Pralits et al. (2017) but is briefly reproduced in Appendix A for the sake of completeness. Thus, the stability of the flow is investigated in the full $(\alpha, \beta, Re)$ domain, for $\lambda _{\parallel } = 0.03$ and $\theta =45^\circ$. For a parametric study over the influence of $\lambda _{\parallel }$ and $\theta$ on the linear stability analysis, the reader is referred to Pralits et al. (2017). Mainly, two main unstable zones can be identified: a patch for small values of $\alpha$ and $\beta \ne 0$ and a horseshoe region at small $\alpha$ and $\beta$ closer to zero. Modes with a large negative spanwise wavenumber have significantly larger growth rates. Note that this unstable region was also observed in the work of Zhai et al. (2023). A strong asymmetry of the neutral isosurface is observed: perturbations are most unstable in the opposite direction of the cross-flow base flow velocity ($\beta < 0$). This has already been noted for swept wings (Mack 1984). Slices of the neutral isosurface in the $\alpha \unicode{x2013}Re$ and $\alpha \unicode{x2013}\beta$ planes, showing the regions of strongest instability, are provided in figure 2.

Figure 2. Neutral curve in the $\alpha \unicode{x2013}\beta$ plane for $Re=12\,000$ (a) and $\alpha \unicode{x2013}Re$ plane for $\beta =-0.5$ (b). Black contours correspond to the projection of the neutral isosurface $\omega _i(\alpha,\beta, Re)=0$ in the respective planes.

For further insight, figure 3 displays representations of the eigenfunctions of the most unstable mode for both regions. The horseshoe instability region is reminiscent of the unstable region of both classic and homogeneous slip Poiseuille flows. The usual shape of a 3-D TS wave (Zang & Krist 1989) is also retrieved. A similar result was found in the analyses of Chai & Song (2019) and Xiong & Tao (2020) in which the linear stability of a channel flow with anisotropic slip boundary conditions was investigated. The boundary conditions that was considered in their works are equivalent, in the present formalism, to setting $\theta =0^\circ$ or $\theta =90^\circ$, yielding a non-zero spanwise slip but a zero spanwise component for the base flow. Thus, the three-dimensionality of the TS wave seems to stem from the presence of spanwise slip and is not a cross-flow-related effect.

Figure 3. Eigenfunctions of the most unstable mode for the two identified regions of interest: for (a,b), $\alpha = 0.7$, $\beta =-0.6$ and $\omega = 0.0004-0.1643i$ (3-D TS wave) and for (c,d), $\alpha = 0.2$, $\beta =-6$ and $\omega = 0.0033+0.02245i$ (cross-flow vortices). (a,c) Absolute value of the disturbance velocity components. (b,d) The contour plot represents the streamwise velocity disturbance while quiver plot shows the $v$–$w$ cross-flow.

The second region exhibits modes with tilted counter-rotating vortices. These vortices are quasi-stationary ($\omega _i = 0.022$) and propagate almost perpendicularly to the streamwise direction ($\phi = \arctan (\beta _0/\alpha _0) \approx 88^\circ$). These two characteristics are reminiscent of unstable modes of swept boundary layers (Mack 1984) or rotating disc flows (Lingwood 1995). However, in the present configuration, and unlike swept flows, the spanwise velocity component of the base flow does not present an inflection point. This is quite surprising since cross-flow instability is usually depicted as a mainly inviscid mechanism (Saric, Reed & White 2003), thus requiring an inflection point in the spanwise velocity $W_0$. A discussion of the possible role of the spanwise component of the base flow on this instability is reported in Appendix D.

The competition between these two families of unstable modes may depend on the configuration (slip length $\lambda _{\parallel }$ and texture orientation $\theta$) considered: moderate and large slip lengths decrease the wall-normal velocity gradients, thus hampering the development of the viscosity-induced TS waves. In turn, the angle $\theta$ affects the slip direction and the magnitude of the spanwise velocity component, inducing cross-flow related effects, leading to the destabilisation of a mode with $\alpha \approx -\beta$. Note that in the present work the angle $\theta$ has been chosen to maximise the magnitude of the spanwise velocity component and, thus, the cross-flow related effects.

3.2. Non-modal stability analysis

In subcritical conditions, non-modal mechanisms, linked to the non-normality of the NS equations, can induce a transient growth of the energy of small perturbations. Such a growth can be several orders of magnitude higher than the initial energy perturbation (Gustavsson 1991; Reddy & Henningson 1993) thus playing a crucial role in the transition to turbulence. This behaviour is typically investigated by maximising the finite-time amplification of an initial velocity perturbation $\boldsymbol {u}_0$ (Butler & Farrell 1992). Mathematically, the quantity to be optimised can be written as

(3.5)\begin{equation} G(T) = \max_{\boldsymbol{u}_0} \frac{E(\boldsymbol{u}(T))}{E(\boldsymbol{u}_0)} \end{equation}

where $E(T)$ is the kinetic energy of a perturbation at target time $T$, defined as

(3.6)\begin{equation} E(\boldsymbol{u}(t)) = \frac{1}{2L_xL_z} \int_V |\boldsymbol{u}|^2 \,{\rm d}V, \end{equation}

where $V$ is the volume of the computational domain of streamwise and spanwise length $L_x$ and $L_z$, respectively. Following Schmid & Brandt (2014), transient growth is obtained through a singular value decomposition (SVD) of the linearised NS operator:

(3.7)\begin{equation} G(t) = || \exp(\boldsymbol{J}t) ||_E^2 = ||\boldsymbol{V}\exp{(\boldsymbol{\varLambda}t)}\boldsymbol{V}^{{-}1}||_E^2, \end{equation}

where $\boldsymbol {J}$ is the Jacobian of the linearised system and $(\boldsymbol {\varLambda }, \boldsymbol {V})$ are eigenmodes found with the previous linear stability analysis. The code is validated in the case of a smooth channel flow (Schmid & Henningson 2001). With $n=200$ grid points in the wall-normal direction, optimal gain and time of the previous reference are recovered with an error of less than $0.2\,\%$.

Figure 4(a) provides the optimal energy gain $G$ in the $\alpha \unicode{x2013}\beta$ plane, which appears to be strongly affected by the presence of the wall slip. The maximum gain can be obtained for a small but non-zero $\alpha$ and for $\beta =-2$, differently from the channel flow with no-slip walls (Butler & Farrell 1992). Amplification curves corresponding to a maximum energy growth in the smooth channel case are retrieved and compared with those found with anisotropic slip in figure 4(b). The slip setting appears to slightly increase the maximum amplification. No effect on the time at which the optimal gain is reached, $T_{opt}$, can be observed. Note that, for $Re=3000$, the flow in the slip configuration becomes linearly unstable. Previous studies (Min & Kim 2005; Picella et al. 2019) have assessed that isotropic slip conditions have a marginal effect on transient growth. Thus, we conjecture that transient growth is mainly affected by cross-flow related effects. Note that the region of maximum energy growth overlaps with the cross-flow instability region, hinting at the coexistence of both mechanisms. Similar results were found in the non-modal analysis conducted by Breuer & Kuraishi (1994) and Corbett & Bottaro (2001) on swept wings.

Figure 4. (a) Optimal energy gain contours in the $\alpha \unicode{x2013}\beta$ plane for $Re=1000$. The Reynolds number is chosen sufficiently low to ensure stability of the flow. The black dot corresponds to the maximum energy gain obtained for $\alpha =0.1$ and $\beta = -2.0$. (b) Amplification curves for $(\alpha,\beta ) = (0.1,-2)$ in the slip configuration ($\circ$) and for $(\alpha,\beta ) = (0,2)$ in a smooth channel flow ($\square$) for several Reynolds numbers $Re$.

The optimal perturbation ensuring the maximum transient growth is given in figure 5. The initial perturbation (left frame) consists of quasi-streamwise tilted counter-rotating vortices. The vortices are tilted in the opposite direction of the cross-flow component of the base flow and distorted by the slip conditions in the near-wall regions. The cross-flow components have an amplitude one order of magnitude higher than the streamwise component. At the optimal time, the perturbation takes the form of velocity streaks slightly tilted in the cross-flow direction. Non-modal growth is usually the product of two mechanisms: the Orr mechanism and the lift-up effect (Brandt 2014). Since the Orr mechanism is more efficient for short (${O}(10)$) time scales and high streamwise wavenumbers (Butler & Farrell 1992), in the present case the transient growth is most probably linked to the lift-up effect (Landahl 1980). However, the mechanism is modified by the presence of the spanwise velocity component of the base flow. This influence is made clear by following the work of Ellingsen & Palm (1975) who derived a model which showed algebraic growth for the amplitude of streamwise independent perturbations in time. A generalisation of this model for 3-D base flows is not trivial, except for small angles for which the mechanisms remain similar. However, in the specific case of a constant spanwise component of the base flow, an extension of this model is possible and has been derived in Appendix B. According to this model, the streamwise perturbation amplitude can be expressed as follows:

(3.8)\begin{equation} u(t) = u_0 \cos{(\beta W_0t)} - v_0U't. \end{equation}

Figure 5. Optimal perturbation at time $t = 0$ (a) and at target time $t=T = 82$ (b). The contour plot represents the streamwise component of the perturbed velocity whereas the quiver plot depicts the $v$–$w$ cross-flow. Eigenfunctions have been normalised by the maximum of the spanwise velocity at the initial time whereas the maximum of streamwise velocity is used at the optimal time.

While the classical algebraic growth is retrieved, a new term, due to cross-flow effects, appears. This term induces, for early times, an oscillation of frequency $\beta W_0$ of the amplitude of the streamwise velocity of the perturbation.

4. Direct numerical simulations

To better understand the mechanisms arising during the laminar–turbulent transition, supercritical transition of the flow in this configuration is investigated. Such a transition can be triggered by superposing on the base flow (2.6a–c) an unstable wave found from the linear stability analysis of the previous section, say

(4.1)\begin{equation} \boldsymbol{U}(\boldsymbol{x},t=0) = \boldsymbol{U}_0(\boldsymbol{x}) + A\boldsymbol{u}_1(\boldsymbol{x}) , \end{equation}

where $\boldsymbol {u}_1$ is an unstable mode found from the previous linear stability analysis and $A$ is its initial amplitude. In the following, the amplitude is set such that the initial energy of the perturbation is equal to $10^{-5}$ for ensuring a linear phase of exponential growth. Two transition scenarios are considered: one initiated with a 3-D TS wave having $(\alpha _0, \beta _0) = (0.7,-0.6)$, presented in § 5, and the other with a cross-flow mode having $(\alpha _0, \beta _0) = (0.2,-6)$, discussed in § 6. This choice of wavenumbers results from a compromise between their distance to the neutral curve and the numerical cost of the simulation. For the cross-flow scenario, the selected wavenumber is also chosen not to be in the region where transient growth is significant in an attempt to discriminate the different instability mechanisms. A temporal DNS framework is considered and thus, for the 3-D TS, the numerical domain is set to fit exactly one spatial wavelength of the considered perturbation in the streamwise and spanwise directions. As a consequence, only waves with wavenumbers of the form $(n\alpha _0, m\beta _0)$ with $n,m \in \mathbb {Z}$ can develop. For the scenario initiated with the cross-flow mode, the large spanwise wavenumber leads to a small computational box in the spanwise direction which was deemed too constraining. Instead, a larger numerical domain is used to avoid restrictions on the wavelengths of developing modes. Ultimately, the domain sizes are $[L_x, L_y, L_z] = [2{\rm \pi} /\alpha _0, 2, 2{\rm \pi} /\beta _0]$ and $[L_x, L_y, L_z] = [2{\rm \pi} /\alpha _0, 2, 6{\rm \pi} /\beta _0]$, respectively. In both cases, the Reynolds number is $Re=12\,000$ and a constant flow rate in the streamwise direction is imposed. The numerical parameters employed in these two cases can be retrieved in table 1.

Table 1. Numerical parameters for the two simulations, where $N_p$ denotes the polynomial order of the elements.

The simulations have been performed with the spectral element incompressible solver Nek5000 (Nek5000 Version 19.0 2019). The code is based on a splitting method which leads to solving a set of time-dependent problems and a final correction step together with a $P_n-P_{n-2}$ spatial discretisation. Convective terms are treated with an explicit extrapolation scheme of order three whereas viscous terms are solved with a backward differentiation scheme (also of order three). In the standard Nek5000 distribution, essential (Dirichlet) and natural (Neumann) boundary conditions are implemented. The homogenised Robin boundary conditions (2.3) must be implemented. The full details of this implementation can be found in the PhD dissertation of Picella (2019). Spectral elements in the wall-normal direction have been distributed following a Chebyshev grid to further increase the spatial resolution near the walls.

For further insight, several indicators are monitored through the simulation such as the spectral energy associated with Fourier modes defined as (Zang & Krist 1989; Schmid & Henningson 2001)

(4.2)\begin{equation} \hat{E}_{k_x, k_z} = \tfrac{1}{2}\int_{{-}1}^1 |\hat{u}_{k_x,k_z}(y,t)|^2\,{{\rm d}y}, \end{equation}

where $\hat {u}_{k_x,k_z}(y,t)$ is the Fourier mode of the perturbation velocity field with streamwise and spanwise wavenumbers $k_x$ and $k_z$, respectively. Following the formalism of Schmid & Henningson (2001), the mode denoted as $(m,n)$ corresponds to the $(m \alpha _0, n \beta _0)$ Fourier mode. Transition to turbulence will also be investigated through the evolution of the kinetic energy density of the disturbance (as defined previously) and the friction Reynolds number, which is defined as

(4.3)\begin{equation} Re_{\tau}= \sqrt{Re|\partial_y \bar{U}(\boldsymbol{x},t)|_{\pm1}}, \end{equation}

where the overbar denotes spatial averaging on the $x$–$z$ plane. This quantity will sharply increase during the turbulent breakdown while remaining almost constant for both laminar and turbulent flows. In the laminar case, the friction Reynolds number can be found analytically and reads

(4.4)\begin{equation} Re_{\tau} = \sqrt{\frac{2Re}{1+3L_{11}}} \approx 150. \end{equation}

In addition to the two supercritical transition scenarios described previously, two other simulations have been performed: a channel flow with no-slip boundary conditions to serve as a reference and a channel flow with anisotropic slip but no cross-flow (corresponding to the $\theta =0^\circ$ case). These two simulations have been performed starting from the same initial energy $E_0 = 10^{-5}$ and with the exact same numerical domain than in the TS scenario. An overview of the different cases for the kinetic energy and friction Reynolds number indicators can be seen in figure 6. As expected, slip significantly delays the transition and reduces friction in the vicinity of the wall. Note that for $\theta =0^\circ$, transition could not be observed, underlining the destabilising effect of the spanwise velocity base flow. These scenarios are discussed in much more detail in the following two sections.

Figure 6. Evolution of the disturbance kinetic energy density (a) and $Re_{\tau }$ (b) as a function of time for different texture parameters and transition scenarios.

5. First scenario: 3-D TS waves

5.1. Overview of the transition

The transition scenario initiated by a 3-D TS wave with $(\alpha _0, \beta _0) = (0.7,-0.6)$ is investigated. The time evolution of the kinetic energy evolution and friction Reynolds number are provided in figure 7, whereas figure 8 depicts the evolution of the spectral energy for several Fourier spatial modes. In the initial phase, the kinetic energy density displays an exponential increase, with a growth rate in perfect agreement with the linear stability analysis. Figure 8 shows that the Fourier mode $(1,1)$, which is the only finite-amplitude mode present in the flow at early times, follows this evolution. Figure 9 (top row) shows that at $t=200$, the perturbation is still characterised by a TS wave shape.

Figure 7. Evolution of the disturbance kinetic energy density (a) and $Re_{\tau }$ (b). The dashed line has a slope $2\omega _r$, where $\omega _r \approx 0.00044$ is the growth rate of the most unstable mode found by means of linear stability analysis.

Figure 8. (a) Time evolution of the Fourier modes $(1,1)$, $(0,16)$ and $(0,1)$. The dashed line has a slope $2\sigma \approx 0.0406$. Note the oscillations in the early evolution of mode $(0,16)$. (b) Time evolution of the Fourier modes ranging from $(0,1)$ to $(0,9)$: blue, $(0,1)$; purple, $(0,2)$; brown, $(0,3)$; light green, $(0,4)$; grey, $(0,5)$; dark green, $(0,6)$; red, $(0,7)$; orange, $(0,8)$; pink, $(0,9)$. All these modes appear to experience a significant energy growth during the second phase of the transition process.

Figure 9. Snapshots of the flow at different times $T$. From left to right and top to bottom: $T=200, 400,600,800,900$. Isosurfaces of the $\lambda _2$-criterion, $\lambda _2=-10^{-5}$ for the two first rows, $-10^{-4}$ for the third and fourth rows and $0.25$ for the last row. The contours depict the streamwise velocity at the wall. The flow is from bottom to top, and left to right. For the sake of clarity, only half the channel is shown.

Meanwhile, at $t \approx 100$, the mode $(0,16)$ starts oscillating with a period $T_0 = 47$. The period of the oscillations is in good agreement with that analytically predicted from (3.8), namely $2{\rm \pi} /(\beta W_0) = 46.58$, hinting at the presence of the modified lift-up described previously. Note that these oscillations can be observed only on this high-wavenumber mode, whereas other streamwise-invariant Fourier modes such as $(0,1)$ (see figure 8) does not display an initial oscillatory phase. A possible reason for this behaviour is that both the target time at which maximum transient growth is recovered, $T_{opt}$, and the oscillation period, $T_0$, are much larger for small spanwise wavenumbers. Furthermore, for $(0,16)$, it can be found that $T_{opt} = 100$ and $T_0 = 47$ whereas $T_{opt} = 615$ and $T_0 = 639$ for the $(0,1)$ mode, which is too large a period to be observed during the linear phase of the perturbation evolution.

After $t\approx 280$, the $(0,16)$ mode exhibits exponential growth, characteristic of a secondary instability. At $t\approx 520$, the $(0,16)$ mode becomes more energetic than the primary one. At the same time, a sharp increase in the kinetic energy indicates that the first step of transition to turbulence has initiated. This translates into a change in the topology of the mode characterised by the onset of near-wall streamwise-elongated coherent structures, shown in figure 9 (second and third rows from the top). These streaky structures are characterised by a large spanwise wavenumber, $\beta =-9.6$, corresponding almost exactly to $16$ times the fundamental spanwise wavenumber $\beta _0$ (namely, $\beta =16\beta _0= -9.6$).

In this transition phase, energy is taken from the 3-D TS wave and transferred to these oblique waves (see figure 9, second and third rows from the top) ultimately leading to the disappearance of the TS wave. As can be seen in figure 10, showing the spatiotemporal evolution of the streamwise (left) and spanwise (right) components of the velocity, these streaky structures are not streamwise independent and are oriented with an angle $\phi = 9^\circ$ with respect to the streamwise direction, almost perpendicular to the initial 3-D TS wave. Note also how the first phase of the transition leads to the decrease of the spanwise velocity component, confirming the streaky nature of the structures.

Figure 10. Spatiotemporal evolution of the streamwise (a,c) and spanwise (b,d) components of the velocity at $x=0$ and $y=-0.85$. For visualisation purposes, the first (second) phase of the transition is depicted separately in the first (second) row.

At $t = 650$, the kinetic energy reaches a plateau for $E = 3\times 10^{-3}$ as the secondary instability saturates (see figure 7). Ultimately, tertiary instability triggers, at $t\approx 800$, a dramatic increase in the friction Reynolds number, indicating the breakdown into a fully turbulent state. This final transition can be linked to the $(0,1)$ mode, as suggested by its sudden growth at such time. Note that this growth is not specific to the $(0,1)$ mode and seems to be shared by all the streaky modes from $(0,1)$ to $(0,9)$ (see figure 8), demonstrating that this behaviour is somewhat general for the streamwise-independent modes and not strictly linked to the domain size. The nature of this final instability is quite complex: it is highly nonlinear as indicated by the growth of all these spatial Fourier modes. It also and appears to be, to some extent, related to streak instability. Figure 10(c,d) depicts streaky structures which, quite interestingly, have lost their orientation and are now aligned with the streamwise direction. The streak instability can be observed at time $t\approx 700$ in the two bottom rows of figure 9 and in the bottom left frame of figure 10.

Quite notably, this transition scenario is fundamentally different from the usual transition scenario in channel flows. Indeed, supercritical transition to turbulence in channel flow is well documented in the no-slip case (Zang & Krist 1989; Schmid & Henningson 2001): after an initial exponential growth phase, secondary instability of a subharmonic mode appears. The interaction between the TS wave and the secondary modes leads to a peak–valley structure (Herbert 1985; Asai & Nishioka 1989), which, in turn, results in the formation of a staggered pattern of $\lambda$-vortices. Then, $\lambda$-vortices subsequently develop into hairpin vortices that will ultimately breakdown into turbulence. In the slip case, with a homogeneous Navier slip boundary condition, Picella et al. (2019) have shown that a delay of the laminar–turbulent transition is possible, but the nature of the instability mechanisms remains virtually the same.

In the present case, the shear misalignment of the base flow appears to have a dramatic effect on the transition scenario. Similarly to what we have reported concerning the primary instability, it might originate a new secondary instability mechanism, which would be investigated in the next subsection. If these cross-flow-related secondary modes have higher growth rates than the subharmonic modes leading to $\lambda _2$ vortices (Herbert 1985), the secondary phase of transition would be dominated by cross-flow-related mechanisms.

As observed in swept flows (Serpieri 2018), during this secondary phase, oblique waves appear and eventually, through nonlinearities, saturate. The saturated cross-flow vortices contain several strong shear layers, which can, under certain circumstances, further destabilise. The most common instability mechanism observed in such flows (Bippes 1990, 1999) involves the shear layer at the bottom of the vortex, created from the circulation of high-speed fluid towards lower velocities region near the wall.

In the present case, the slip boundary conditions reduce the efficiency of this mechanism by lowering the wall-normal velocity gradients near the walls. Instead, the flow keeps accelerating in the spanwise direction until it reaches a neighbouring vortex, thus creating stagnation points and strong shear layers in the spanwise direction. Transient growth in this region leads to the formation of streaks (Guegan, Huerre & Schmid 2007) aligned with the streamwise direction. Streak instability concludes the transition to turbulence.

This scenario of transition will be corroborated in the next subsection, by means of secondary stability analysis and comparison with DNS.

5.2. Floquet stability analysis

In an effort to provide a more qualitative description of the secondary phase of the transition, a secondary stability analysis, based on Floquet theory, is carried out. This method has been applied to both channel flows (Herbert 1983) and 2-D boundary layers (Herbert 1985) for which it has successfully predicted the secondary instability of a 2-D TS wave. An extension to 3-D base flows was proposed and applied most notably to swept wings (Fischer & Dallmann 1991; Janke & Balakumar 2000; Liu, Zaki & Durbin 2008). Secondary stability equations in a velocity-vorticity formulation can be found in Fischer & Dallmann (1991). A primitive variables formulation for 2-D base flows only, can be found in Schmid & Henningson (2001).

In this work, due to the three-dimensionality of the problem under consideration, we use a formulation of Floquet theory and subsequent modal expansion somewhat different (but nevertheless equivalent) from those used in these articles. In the study of Herbert (1983), the Squire theorem ensures a 2-D initial perturbation while, for swept flows cases, an initial 3-D perturbation should be considered. For this reason, Floquet theory for swept flows is usually applied in a rotated frame aligned with the 3-D wave direction of propagation, thus effectively reducing the base flow to two dimensions. Precisely, the base flow for secondary stability analysis is constructed as a superposition of an unstable mode and the laminar flow profile such as

(5.1)\begin{equation} \boldsymbol{U}_1(x,y,t) = \boldsymbol{U}_0(y) + A\boldsymbol{\hat{u}}^{TS}(y)\exp({{\rm i}\alpha_0x + {\rm i}\beta_0 z -{\rm i}\omega_0t}) \end{equation}

with $A$ the amplitude of the unstable 3-D TS wave obtained from linear stability analysis. Considering that the TS wave moves at the velocity $c_x$ and $c_z$ in the $x'$ and $z'$ directions, respectively, the time dependence of the base flow is removed when the moving frame $(x',z') = (x-c_xt, z-c_zt)$ is considered. Linearising the NS equations around this new base flow in the co-moving frame leads to a system of ordinary differential equations (see Appendix C) with $2{\rm \pi} /\alpha _0$-periodic and $2{\rm \pi} /\beta _0$-periodic coefficients in the $x$- and $z$-directions, respectively. Thus, the Floquet theorem is applied in both the streamwise and spanwise directions:

(5.2)\begin{equation} \boldsymbol{q_1}(x',y,z',t) = \boldsymbol{\tilde q}(x',y,z') \operatorname{e}^{\gamma x'} \operatorname{e}^{\mu z'} \operatorname{e}^{\sigma t} \end{equation}

with $\gamma = \gamma _r + {\rm i}\gamma _i$ and $\mu = \mu _r + {\rm i}\mu _i$ the Floquet parameters in the streamwise and spanwise directions. Then, $\boldsymbol {\tilde q}(x',y,z')$ is Fourier transformed in both the $x'$- and $z'$-directions and introducing the detuning factors $\epsilon = \gamma _i/\alpha _0$, $\delta = \mu _i/\beta _0$, the general form of the solution can be expressed as

(5.3)\begin{align} \boldsymbol{q_1}(x',y,z',t)& = \exp({\sigma t})\exp({\gamma _r x'}) \exp({\mu_r z'}) \sum_{m,n={-}\infty}^{+\infty} \boldsymbol{\tilde q}_{m,n}(y) \exp({\rm i}\alpha_0(m + \epsilon)x' \nonumber\\ &\quad + {\rm i}\beta_0(n + \delta)z') . \end{align}

When a numerical solution is sought, the modal expansion (5.3) needs to be truncated, usually with the lowest possible number of modes. Unfortunately, reaching the spanwise wavenumbers observed in the DNS would require a large number of modes in the $z'$-direction, which make the problem too computationally expensive. In order to relax this constraint, the modal expansion in the spanwise direction is reduced to only one mode by fixing $n$ and $\delta$ and introducing $\beta _1 = \beta _0(n + \delta )$, which represents the effective secondary spanwise wavenumber. Since $n \in \mathbb {Z}$ and $-1/2 < \delta \leq 1/2$, the quantity $n+\delta$ spans all real numbers bijectively such that every spanwise wavelength is accessible through a unique choice of $n$ and $\delta$.

Temporal stability is investigated, thus $\gamma _r=\mu _r=0$ and the real part of $\sigma$ indicates the growth rate of the secondary instability. Finally, the modal decomposition (5.3) reduces to

(5.4)\begin{equation} \boldsymbol{q_1}(x',y,z',t) = \exp({\sigma t})\exp({{\rm i}\beta_1 z'})\sum_{m={-}\infty}^{+\infty} \boldsymbol{\tilde q}_{m}(y) \exp({{\rm i}\alpha_0(m + \epsilon)x'}) . \end{equation}

Introducing (5.4) into the linearised NS equations leads to an infinite set of equations which, once truncated, can be recast in an eigenvalue problem likewise the primary stability problem.

All the results given here are obtained with the lowest possible truncation that is ${m=0,1}$ in the subharmonic ($\epsilon =1/2$) cases and $m=-1,0,1$ for the fundamental ($\epsilon = 0$) and detuned ($0 < \epsilon < 1/2$) modes. For the sake of clarity, the derivation of the secondary stability equations and the eigenvalue problem is fully detailed in Appendix C.

In both classic channel and boundary layer flows, fundamental and detuned waves play a limited role in the transition process (Herbert 1983, 1985) as subharmonic instabilities have stronger growth rates. However, this might not be true in our configuration, due to the strong asymmetry present in the base flow.

The dependence of the growth rate of the most unstable mode on the detuning factor is investigated in figure 11, for $Re=12\,000$. A strong secondary growth rate can be reached for detuned modes and large spanwise wavenumbers, in the range of those recovered for the primary cross-flow instability. The maximum secondary growth rate is $\sigma _r = 0.0230$ for $\epsilon = 0.23$ and $\beta _1=-8$ but the instability region is quite large and similar growth rates can be reached for larger $\beta _1$. Precisely, the growth rate corresponding to the spanwise wavelength observed in the DNS (namely, $\beta = -9.6$), is equal to $\sigma _r = 0.020$.

Figure 11. (a) Growth rate $\sigma _r$ of the most unstable mode as a function of the spanwise wavenumber $\beta _1$ for ${\epsilon =0}$ (fundamental), $\epsilon = 0.22$ (detuned) and $\epsilon = 1/2$ (subharmonic). (b) Contour plot of the growth rate $\sigma _r$ of the most unstable mode for $Re=12\,000$ as a function of both the detuning factor $\epsilon$ and the spanwise wavenumber $\beta _1$. Black dot corresponds to the maximum growth rate obtained for $\epsilon =0.23$ and $\alpha =-8$ whereas the black line denotes the $\sigma _r = 0.02$ isocontour (corresponding to the growth rate found from the DNS). The broken line denotes the $\sigma _r = 0.015$ isocontour and allows the identification of secondary maxima.

The spectrum obtained for $Re=12\,000$, $\beta _1 = -9.6$, $\epsilon = 0.23$ and $A = 10^{-5}$, shown in figure 12, exhibits three unstable and one marginally unstable mode. All the modes share a similar structure: the $m=0$ mode, yielding an effective streamwise wavenumber $\alpha = \epsilon \alpha _0 = 0.16$, is dominant in comparison with the $m=-1$ and $m=1$ modes. The wavenumber in the streamwise direction is too small to be retrieved in the computational domain thus explaining, to some extent, the onset of streaky structures. The eigenvalues associated to these modes are summarised in table 2. Considering the dominant $m=0$ mode, frequencies in the laboratory frame are retrieved through the following equation:

(5.5)\begin{equation} \sigma_i^{0} = \sigma_i^{TS} - \beta_1 c_z - \alpha_0 \epsilon c_x. \end{equation}

Figure 12. Secondary instability spectra for $Re=12\,000$, $\beta _1=-9.6$, $A = 10^{-5}$ and $\epsilon =0.23$. Eigenvalues are coloured according to their respective Fourier mode. The values of the frequency and growth rate of the four most unstable modes are given in table 2.

Table 2. Most unstable eigenvalues of the secondary stability analysis for $Re=12\,000$, $\beta _1 = -9.6$, $\epsilon = 0.23$ and $A = 10^{-5}$. For a single mode, two frequencies, $\sigma _i^{TS}$ and $\sigma _i^{0}$, are given: the former is the frequency of the mode in the frame moving with the TS wave whereas the latter corresponds to the frequency in the laboratory frame.

For the four most unstable modes, secondary perturbations are reconstructed based on (5.3). A slice of these modes in the $y\unicode{x2013}z$ plane is provided in figure 13, showing that they share a rather similar structure. In the spanwise direction, one can see alternated high and low streamwise-velocity patches. The streamwise velocity component is one order of magnitude higher than the cross-flow components, and is strongly localised in the region where the amplitude of the TS wave reaches a maximum. Cross-flow components form vortical structures concentrated near the walls and similar to those found for cross-flow instabilities. More precisely, mode 1 consists of oblique vortices together with oblique streaky structures for the streamwise velocity component. Modes 2 and 4 are almost identical besides a phase shift. One can also observe a stagnation point near the wall at $z=0.58$ for mode 2. Mode 3 appears peculiar as its vortices sit on top of a region of low velocity. For all these modes, the vortices push high-momentum fluid towards the wall where, due to the slip boundary conditions, it is strongly accelerated in the spanwise direction. When this high-momentum flow reaches the neighbouring vortices, it is re-ejected upwards back in the flow, creating stagnation points and strong spanwise shear layers in the process.

Figure 13. Slice at $x=L_x/2$ of the reconstructed secondary perturbations of the four most unstable modes for $Re=12\,000$, $\beta _1 = -9.6$, $\epsilon = 0.23$ and $A = 10^{-5}$. From top to bottom and left to right, modes are ordered from the most unstable to the least unstable. The contour plot represents the streamwise velocity disturbance whereas the quiver plot shows the $v$–$w$ cross-flow.

The main features of the most unstable modes recovered by Floquet analysis are compared with those observed in the DNS in table 3. For the DNS, the secondary growth rate is extracted from figure 8. A good agreement can be found between the values. A comparison of the spatial structure of the secondary perturbations obtained from the DNS and some of the modes obtained by secondary stability analysis is provided in figure 14. Snapshots of the flow are taken at $x=L_x/2$ and $t=500$. Both the base flow and the initial TS wave are subtracted from the snapshot to retrieve the secondary perturbations. The snapshot is then divided into portions of size $2{\rm \pi} /\beta _1$ to isolate vortical structures.

Table 3. Comparison between spanwise wavenumbers and growth rates obtained from the DNS and from secondary stability analysis.

Figure 14. Comparison between snapshots extracted from the DNS (a,c) and secondary stability modes (b,d): mode 1 (a,b) and mode 3 (c,d). Secondary perturbations from the DNS are extracted from snapshots taken at $t=500$, by subtracting the base flow and the primary disturbance. The contour plot represents the total streamwise velocity disturbance whereas the quiver plot shows the $v$–$w$ cross-flow.

The spatial structures compare rather well, although for mode 1 (top row), one can observe that in the DNS the counter-rotating vortices are located rather farther from the wall. Regarding mode 3, its location in the wall-normal direction can be retrieved despite the size of the vortex is slightly smaller. These discrepancies could be explained from different factors: as can be seen, secondary stability analysis yields a range of unstable spanwise wavenumbers which includes the main wavenumber observed in the DNS. The growth rates are in perfect agreement. The small discrepancy found between the frequencies could arise from the fact that the shape assumption is likely not completely valid and the mean flow is slightly distorted during the transition. In addition, only the frequency of the dominant mode $m=0$ was considered whereas, formally, each mode of the Floquet expansion has its own frequency in the laboratory frame.

6. Second scenario: cross-flow modes

6.1. Overview of the transition

The second transition scenario is initiated injecting on the base flow the cross-flow unstable mode with $(\alpha _0, \beta _0) = {(0.2,-6)}$. The time evolution of the kinetic energy evolution and friction Reynolds number is provided in figure 15, whereas figure 16 depicts the evolution of the spectral energy for several Fourier spatial modes along with spatiotemporal plots.

Figure 15. Evolution of the disturbance kinetic energy density (a) and $Re_{\tau }$ (b) for the cross-flow transition scenario. The dashed line has a slope $2\omega _r$ where $\omega _r \approx 0.003267$ is the growth rate of the most unstable perturbation found with the linear stability analysis.

Figure 16. (a) Time evolution of the Fourier modes $(1,1)$, $(2,0)$, $(7,0)$ and $(0,1)$ for the cross-flow transition scenario. The dashed line corresponds to the exponential primary growth rate. (b) Spatiotemporal evolution of the streamwise component of the velocity at $x=0$ and $y=-0.85$.

At first, exponential growth of the disturbance kinetic energy is observed in figure 15. Note that the slope of the curve corresponds to twice the growth rate $\omega _r = 0.0033$ retrieved from linear stability analysis. After $t=600$, the kinetic energy starts saturating at $E \approx 1.8\times 10^{-3}$. This saturation most likely arises from stabilising nonlinear effects and is also observed in a cross-flow-induced transition for swept flows. As shown in the first snapshot of figure 17, apart from slight deformations near the walls due to the boundary conditions, the shape of the quasi-stationary ($\omega _i = 0.0225$) cross-flow vortices does not evolve much during this phase.

Figure 17. Snapshots of the flow at $T=900,1010,1030,1060,1080,1150$ (from top to bottom and left to right). Isosurfaces of the $\lambda _2$-criterion, $\lambda _2=-10^{-5}$ (a,b), $\lambda _2=-10^{-4}$ (c,d) and $\lambda _2=-10^{-1}$ (e, f), respectively, and contours of the streamwise velocity at the wall. For the sake of clarity, only half the channel is shown.

In the meantime, it can be seen from figure 16(a) that the mode $(2,0)$ also experiences strong energy growth. Since the mode is linearly stable and is not expected to undergo strong transient growth, only nonlinear mechanisms could explain this growth. Nonlinear forcing of the mode $(2,0)$ by the initial perturbation is possible as it can be observed in an oblique waves transition scenario for a boundary layer (Schmid & Brandt 2014). Usually, the $(2,0)$ mode is damped and not very receptive to forcing, meaning it is not relevant in the transition scenario. For 3-D flows, El-Hady (1989) demonstrated that resonant wave triads play an important role in the stability due to the large number of interaction between the possible instabilities.

At $t = 1000$, the top of the vortices start oscillating in the streamwise direction (see figure 17b,c). This can also be seen on the spatiotemporal plot of figure 16(b), where one can note the much higher frequency of the secondary instability. This secondary instability can be traced back to the Fourier mode $(7,0)$ in figure 16(a). Secondary growth rate can be extracted from the slope of its time evolution and is equal to $\sigma _r \approx 0.084$.

At $t \approx 1050$ (see figure 17d,e), spanwise-independent coherent structures appear as shown in figure 18. These structures appear on top of the cross-flow vortices and seem to be tilted. The bottom snapshot of figure 18 is particularly interesting: on top of the wave-like pattern of the initial disturbance, secondary vortices can be observed. These also appear in the near-wall region between two cross-flow vortices. These vortices are most likely a consequence of the slip boundary conditions: as the fluid is accelerated in the spanwise direction near the walls, strong spanwise shear layers, which are susceptible to further destabilise, are created. Energy is transferred from the cross-flow vortices to the secondary perturbation. Distortion of the base flow, in both streamwise and spanwise directions, is also highly likely.

Figure 18. Cross-flow dynamics at different time $T$. From top to bottom: $T=900,1000,1050, 1060,1080$. The cross-section is taken at the streamwise position $x_s = 15.71$ and only $y \in [-1,-0.7]$ is represented. Isocontours of the $\lambda _2$-criterion [$-1\times 10^{-4} \longleftrightarrow -1\times 10^{-2}$]. Contour plot of the $u$-component of the disturbance velocity. Quiver plot depicts the $(v$–$w)$ cross-flow of the perturbation. White contours represent the magnitude of the velocity $[-5\times 10^{-2}, -3\times 10^{-1}; -5\times 10^{-2}]$.

The remaining part of the laminar–turbulent transition is rather complex as several instability mechanisms get intertwined. In the last two snapshots of figure 17, streaks can be seen developing in the near-wall region, causing the energy growth of the $(0,1)$ Fourier mode in the left frame of figure 16. In figure 18, the vortices described previously have combined to form vortex quadrupoles which create strong recirculation in the flow. Ultimately, the breakdown to turbulence seems to be related to the displacement of high-velocity fluid in the upper part of the channel (fourth snapshot of figure 18) towards lower-velocity regions near the wall (last snapshot of figure 18). In the vicinity of these regions turbulent wedges originate (see figure 16b) which quickly propagate to the whole channel.

6.2. Secondary stability analysis of cross-flow vortices

Due to the presence of several strong shear layers in the primary cross-flow vortices, these are highly likely to destabilise. The observation of the transition scenario also strongly hints at the presence of modal secondary instability mechanisms. In the case of swept flows, numerous efforts on the secondary instability of cross-flow vortices have been performed, albeit, to the best of the authors’ knowledge, always for boundary flows with no-slip boundary conditions. While a short introduction to secondary instability of cross flow vortices is produced in the following, the reader is referred to Saric et al. (2003) for a full review.

One of the first key contribution to this issue has been the work of Malik et al. (1999) which classified the secondary modes into two main families. Precisely, type-I unstable modes are linked to the velocity gradient of the streamwise component of the velocity in the spanwise direction, whereas type-II modes originate from gradients in the wall-normal direction. Type-I modes are located in the outer part of the primary vortices whereas type-II are situated on top of them. Type-II also tend to have a higher frequency than type-I modes. A third type of mode was identified in studies by Fischer & Dallmann (1991) and Janke & Balakumar (2000) through Floquet analysis of the cross-flow vortices. This low-frequency type-III mode is linked to nonlinear interactions between stationary and travelling primary modes. These three families of modes have been also experimentally retrieved by Kawakami, Kohama & Okutsu (1999) and White & Saric (2005). Later, using DNS, Wassermann & Kloker (2002, 2003) showed that co-rotating helicoidal structures superimposed on the upwelling region of the primary vortices were characteristics of the type-I mode. These structures would be convected downwards if the unsteady forcing was switched off, thus confirming the convective nature of the secondary instability as already suggested by Kawakami et al. (1999) and Koch (2002).

Bonfigli & Kloker (2007) also studied the development of the secondary instabilities of cross-flow vortices through the combined use of a spatial DNS and secondary linear stability theory (SLST) based on Floquet theory applied on a saturated primary flow (Malik, Li & Chang 1994; Janke & Balakumar 2000; Koch et al. 2000; Koch 2002). Growth rates and structures of the type-I and type-III modes were retrieved and good agreement was shown between the two techniques. SLST predicted the instability of type-II modes but these were not observed in the DNS. In the particular case of travelling vortices, type-III modes were not reported as unstable, which was expected since this mode arises from the generation of spanwise modulation induced by the stationary primary vortices. Another key conclusion of this study is the nature of the mechanism behind the secondary instability: both type-I and type-II are found to be related to Kelvin–Helmoltz instabilities. In all cases, the primary flow on which SLST is carried out is paramount and its acquisition is not straightforward as nonlinearities lead to an important distortion of the base flow. The shape assumption no longer holds and Herbert's secondary stability theory cannot be applied directly.

The primary flow can be retrieved by using parabolised stability equations as in Malik et al. (1999) and Janke & Balakumar (2000), or extracted directly from a saturated DNS (Messing & Kloker 2004; Bonfigli & Kloker 2007). Recently, Groot et al. (2018) based their secondary stability analysis on an experimentally acquired primary flow. In any case, the extraction process is quite complex and several issues arise. In swept boundary layers, the orientation of the vortex axis of the cross-flow vortices evolves with the streamwise direction. Consequently, secondary linear stability depends on the position $x=x_{SLST}$ where the primary flow is extracted. This is not the case in our configuration, as a consequence of the choice of a temporal DNS.

As in Bonfigli & Kloker (2007) and several others, we introduce a vortex-oriented system $(x_v,y,z_v)$ aligned with the vortex axis of the primary disturbance and moving with it. In this frame, the primary flow $\boldsymbol {Q}_v^0(y,z_v)$ is streamwise independent and periodic in the spanwise direction. In the framework of Floquet theory, it is Fourier expanded in the spanwise direction and takes the following modal form:

(6.1)\begin{equation} \boldsymbol{Q}_v^0(y,z_v) = \sum_{n={-}\infty}^{+\infty} \boldsymbol{\hat{q}}^0_{v,n}(y)\exp({{\rm i}n\beta_vz_v}) \end{equation}

with $\beta _v =\sqrt {\alpha _0^2 + \beta _0^2} \approx \beta _0$ the effective wavelength in the spanwise direction $z_v$. In this framework, the continuity equation reduces to

(6.2)\begin{equation} \partial_y v^0_{v,n} + {\rm i}n\beta_vw^0_{v,n} = 0, \quad \forall n, \end{equation}

directly relating the two cross-flow components of the velocity. Thus, a choice must be made between enforcing continuity and extracting both components from the DNS. Ultimately, this might have a non-negligible effect on the secondary growth rates found from the SLST (Malik et al. 1994; Bonfigli & Kloker 2007).

To overcome this and other issues of Floquet analysis, here we use a local 2-D stability (Tatsumi & Yoshimura 1990) approach. Several reasons motivate this choice: first and foremost, 2-D local stability theory encompasses Floquet stability theory. In addition, as shown previously, Fourier decomposition in the spanwise direction of the primary flow leads to some difficulties with the continuity equation. In a 2-D framework, this decomposition is not necessary and these issues do not arise. Ultimately, 2-D stability analysis allows to take into account also instabilities spanning multiple cross-flow vortices, whereas other studies only considered a single pattern of the primary flow for their secondary stability analyses.

For carrying out the 2-D secondary stability analysis, the flow is decomposed into the primary state $\boldsymbol {Q}_v^0(y,z_v)$ and a secondary perturbation $\boldsymbol {q}_1(x_v,y,z_v,t)$. Furthermore, the secondary perturbation can be expanded in the following way:

(6.3)\begin{equation} \boldsymbol{q}_1(x_v,y,z_v,t) = \tilde{\boldsymbol{q}}(y,z_v)\exp({{\rm i}\alpha_v x_v - \sigma t}), \end{equation}

where $\alpha _v$ and $\sigma$ are complex numbers representing the streamwise wavelength and the pulsation in the vortex-oriented reference frame, respectively. Both temporal (Malik et al. 1999; Koch et al. 2000; Wassermann & Kloker 2002) and spatial (Janke & Balakumar 2000; Bonfigli & Kloker 2007) approaches have been considered previously. Koch et al. (2000) found a set of equations relating both frameworks through a generalisation of the Gaster transformation to 3-D flows. For channel flows, a natural choice is to consider a temporal framework. As previously, NS are linearised around this new base flow. Stability equations for 2-D problems can be found, for example, in Loiseau (2014).

The primary flow in the vortex-oriented reference frame is shown in figure 19. It was extracted from the DNS at time $t=900$ and $x_{SLST}=15.71$. This flow field can be compared with the flow visualisations from the experiments of Serpieri (2018) (see figure 4.14) and with the primary flow resulting from the DNS of Bonfigli & Kloker (2007) (see figure 7). The water-wave shape of the streamwise velocity component, characteristic of the cross-flow vortices, is clearly visible. However, the spatial extension of the tip of the wave is smaller than that observed for swept boundary flows, likely because the initial spanwise velocity $W_0$ is weaker. The instability mechanism is similar to that described in Serpieri & Kotsonis (2016): cross-flow vortices generate the circulation of low-momentum flow towards high-momentum regions higher up in the channel and conversely. In the present case, the wall-normal shear is maximum near the walls and not in the vortex core due to the slip boundary conditions. The spanwise velocity gradient reaches a maximum in the low-momentum upwelling region on the outer part of the vortex.

Figure 19. Secondary base flow for the secondary stability analysis. The snapshot is extracted from the DNS at time $t \approx 900$, position $x_{SLST} = 15.71$ and is shown in the vortex-oriented reference frame. (a) Contour plot of the velocity component $U$. Quiver plot of the cross-flow perturbation $(v$–$w)$. (b) Growth rate $\sigma _r$ of the secondary most-unstable modes of cross-flow vortices as a function of the effective streamwise wavenumber $\alpha _v$ in the vortex-oriented frame.

Two-dimensional stability analysis of the previously described primary flow is performed. The evolution of the temporal secondary amplification rates in the vortex-oriented frame as a function of the effective streamwise wavelength is shown in figure 19. The maximum growth rate is equal to $\sigma _r = 0.084$ for $\alpha _v=1.4$. Both the growth rate and the effective streamwise wavenumber are in good agreement with those observed in the DNS. A direct comparison between the characteristics of the secondary instability obtained from the SLST and those extracted from the DNS is provided in table 4.

Table 4. Comparison between secondary instability characteristics obtained from the DNS and from the most-unstable mode of the secondary stability analysis.

All the results from the SLST are obtained in the moving vortex-oriented frame. Values of streamwise wavelengths and frequencies in the laboratory frame can be retrieved by expanding (6.3):

(6.4)$$\begin{gather} \sigma_i^{0} = \sigma_{i,v} - \alpha_v c, \end{gather}$$

(6.5)$$\begin{gather}\alpha_1^{0} ={-}\frac{\alpha_v \beta_0}{k}, \end{gather}$$

(6.6)$$\begin{gather}\beta_1^{0} = \frac{\alpha_v \alpha_0}{k}, \end{gather}$$

where $k = \sqrt {\alpha _0^2 + \beta _0^2}$ and $c$ is the velocity of the cross-flow vortices.

For $Re=12\,000$ and $\alpha _v = 1.4$, two families of eigenmodes could be distinguished and are displayed in figure 20. The first, in the left part of the figure, appears strongly localised on the upwelling of the vortex, a representative feature of type-I modes. Among the several unstable modes, no type-II nor type-III structures could be found. This is not unexpected: type-III modes require interactions with unstable stationary cross-flow vortices. Type-II modes stem from the instability of wall-normal velocity gradients but these are sensibly weaker than their counterparts in swept flows. This could be explained by the combined effect of the absence of an inflection point in the spanwise velocity profile $W_0$ and of the slip boundary conditions. Regarding the second mode, it is slightly less unstable ($\sigma _r = 0.0815$ against $\sigma _r = 0.0827$ for the previous mode) and has a smaller frequency ($\sigma _i = 0.4439$ against $\sigma _i = 0.5319$) but presents a strong resemblance with the left snapshot, second row of figure 17 from the DNS (see also figure 18b,c). The corresponding contour plot in figure 20 also displays the existence of a shear layer in the spanwise direction at the top of the cross-flow vortices, right at the location where the spanwise independent structures will develop. Also of interest is the fact that these two perturbations are fundamental: they have the same spatial periodicity than the secondary base flow.

Figure 20. Most unstable secondary modes found for $Re=12\,000$, $\alpha _v=1.4$ (a,b) corresponding respectively to $\sigma = 0.0827 -0.5319 { i}$ (a) and $\sigma = 0.0815 - 0.4439 {i}$ (b) and for $\alpha _v=1$ (c). Slice in the $(y\unicode{x2013}z)$ plane. The contour plots depict the streamwise velocity perturbation whereas the quiver plots represent the cross-flow dynamics.

The most-unstable mode for $\alpha _v=1$ is shown in figure 20(c) and is similar to the first mode for $\alpha _v=1.4$. It is also a type-I mode but, in the present case, a phase shift appears between each cross-flow vortex. The periodicity of the secondary mode is now three times the periodicity of the primary cross-flow vortices, indicating the existence of detuned modes.

7. Conclusions and outlook

In this paper, the influence of slip anisotropy and shear misalignment on the laminar–turbulent transitions of SH surfaces. Following Pralits et al. (2017), anisotropy of the surface is taken into account through the use of a mobility tensor (Kamrin et al. 2010) generalising the Navier slip condition. This set of homogenised boundary conditions leads to the presence of a cross-flow component in the base flow. The presence of an anisotropic slip deeply affects the linear stability of the flow: a new instability region appears for small streamwise and large spanwise wavenumbers. The underlying instability mechanism was found to be similar to the cross-flow instability observed in swept flows. TS waves could also be retrieved but the most unstable modes are 3-D due to the asymmetry of the neutral curve, induced by the cross-flow. The outcome of this competition between these two types of modes depends on the texture parameters (slip length $\lambda _{\parallel }$ and angle $\theta$) considered: an orientation angle close to $45^\circ$ and large slip lengths promote the development of the cross-flow modes. In the present case, the critical Reynolds number is reached for a cross-flow mode with a large negative spanwise wavenumber. Transient growth has been also considered, showing that the optimal perturbation has the form of oblique streaks with non-zero streamwise wavenumber.

DNS have been performed in two distinct flow cases showing different unstable modes. In the TS-wave-initiated transition, instead of the usual subharmonic bifurcation, streaky structures with large spanwise wavenumber have been found, which ultimately destabilise and break down into turbulence, resulting in a two-phase transition. These streaky structures could also be found for anisotropic slip length without the cross-flow component but did not undergo the second phase of transition. The second phase is characterised by a mechanism which promotes the development of small wavenumbers. The second scenario, initiated with a cross-flow-related unstable mode, depicts a highly nonlinear transition, with the flow saturating towards ‘half-mushroom’ structures as seen in Malik et al. (1994), before transitioning to turbulence. In both transition scenarios, secondary stability analysis has been performed. In particular, Floquet analysis has been used for the TS-wave primary instability, whereas a 2-D stability analysis has been used for the cross-flow mode transition. In both cases, the oriented slip and, to a much lesser extent, the cross-flow due to the anisotropy of the roughnesses, strongly change the transition process with respect to that reported in the case of isotropic slip. This underlines the importance of taking into account anisotropy when studying SH surfaces through homogenised boundary conditions. It has been also shown that coherent structures such as the streamwise vortices or the Kelvin–Helmholtz rollers are retrieved by means of secondary stability analysis, and they have a non-negligible influence on the transition scenarios.

Overall, in the DNS, a similar competition between modes than that found using linear stability analysis, is observed: cross-flow modes tend to have the fastest initial energy growth due to a combination of strong transient growth and nonlinear effects. In contrast, the growth of primary and secondary perturbations in the 3-D TS case is much slower. Yet, transition to turbulence in the 3-D TS case occurs more quickly, as the energy in the cross-flow modes tend to saturate for a large period of time preventing a quick breakdown to turbulence. Interestingly, significant delay of the laminar–turbulent transition could be achieved by engineering the textured surface in such a way that it promotes some more favourable transition scenarios. This would still require numerous other DNS for tackling other canonical transition scenarios such as the K-type or O-type.