2. Problem formulation

The model problem involves flow between two parallel plates extending to ±∞ in the x and z directions. Distributed heating and surface grooves spatially modulate this flow. We limit this presentation to modulations that are the spanwise x-coordinate functions. The formulation can be easily extended to modulations being functions of the streamwise z coordinate – its details are not presented due to their length. The conduit is oriented in a way where the gravity g works in the negative y-direction (figure 3). The fluid is assumed to be incompressible with the density ρ variations modelled using the Boussinesq approximation; it has thermal conductivity k, specific heat c, thermal diffusivity $\kappa = k/\rho c$, kinematic viscosity $\nu $, dynamic viscosity $\mu $ and thermal expansion coefficient $\beta $. The non-dimensional equations of motion are

(2.1a–c)\begin{align}\boldsymbol{\nabla }\boldsymbol{\cdot }\bar{V} = 0,\quad\! \frac{{\partial \bar{V}}}{{\partial t}} + (\bar{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{V} ={-}\boldsymbol{\nabla }p + {\nabla ^2}\bar{V} + P{r^{ - 1}}\theta \bar{j},\quad\! \frac{{\partial \theta }}{{\partial t}} + (\bar{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\theta = P{r^{ - 1}}{\nabla ^2}\theta ,\end{align}

where half of the channel height h is used as length scale, velocity vector $\bar{V} = (u\bar{i} + v\bar{j} + w\bar{k})$ is normalized by the viscous velocity scale ${U_v} = v/h$, pressure p is scaled with $\rho U_v^2$, temperature $\theta $ is scaled with $\kappa \nu /(|{g}|\beta {h^3})$, $Pr = v/k$ is the Prandtl number. The relevant boundary conditions are

(2.2a–d)\begin{align}\bar{V}[{y_L}(x),t] = \bar{V}[{y_U}(x),t] = 0,\quad \theta [{y_L}(x),t] = {\theta _L}(x,t),\quad \theta [{y_U}(x),t] = {\theta _U}(x,t),\end{align}

where the subscripts L and U refer to the lower and upper plates, respectively, and ${y_L}(x)$ and ${y_U}(x)$ describe groove geometry.

Figure 3. Schematic diagram of the flow system. Red and blue colours identify hot and cold sections of the walls.

The above system is closed by specifying constraints either in the form of mean pressure gradients in the x and z directions, that is,

(2.3a,b)\begin{align}{\left. {\frac{{\partial p}}{{\partial x}}} \right|_{mean}} = {\wp _x},\quad {\left. {\frac{{\partial p}}{{\partial z}}} \right|_{mean}} = {\wp _z},\end{align}

or in the form of the mean flow rates in the x and z directions, that is,

(2.3c,d)\begin{equation}{Q_x}{|_{mean}} = {\mathrm{\mathbb{Q}}_x},\quad {Q_z}{|_{mean}} = {\mathrm{\mathbb{Q}}_z}.\end{equation}

In the absence of grooves and heating, the velocity and pressure fields of the reference flow have the form

(2.3e–g)\begin{equation}\boldsymbol{V}(x,y,z) = [0,0,w] = [0,\; 0,\; Re(1 - {y^2})],\quad p(x,y,z) ={-} 2\; zRe,\quad {Q_z} = {\textstyle{4 \over 3}}Re,\end{equation}

where the Reynolds number is defined as $Re = {W_{max}}h/\nu = {W_{max}}/{U_v}$ and ${W_{max}}$ denotes the maximum of the z-velocity component. In the computations, reference flow with a specific Re is selected, and its modified state is determined.

Fourier expansions of the form

(2.4a)\begin{gather}{y_L}(x) ={-} 1 + {B_L}{Y_L}(x) ={-} 1 + {B_L}\sum\limits_{n ={-} {N_G}}^{n = {N_G}} {H_L^{(n)}\,{\textrm{e}^{\textrm{i}n\alpha x}}} ,\end{gather}

(2.4b)\begin{gather}{y_U}(x) = 1 + {B_U}{Y_U}(x) = 1 + {B_U}\sum\limits_{n ={-} {N_G}}^{n = {N_G}} {H_U^{(n)}\,{\textrm{e}^{\textrm{i}n(\alpha x + {\varOmega _G})}}}\end{gather}

describe the surface topographies. In the above, ${Y_L}(x)$ and ${Y_U}(x)$ are the shape functions describing plates’ topographies satisfying conditions

(2.4c)\begin{equation}\textrm{max}[{Y_L}(x)] - \textrm{min}[{Y_L}(x)] = 1,\quad \textrm{max}[{Y_U}(x)] - \textrm{min}[{Y_U}(x)] = 1,\end{equation}

${B_L}$ and ${B_U}$ are the peak-to-bottom groove amplitudes at the lower and upper plates, respectively, $H_L^{(n)}$ and $H_U^{(n)}$ are the coefficients of the Fourier expansions describing the shape of grooves at the lower and upper plates, respectively, ${N_G}$ is the number of Fourier modes required to describe the geometry, $\alpha $ is the modulation wavenumber and ${\varOmega _G}$ stands for the phase shift between the upper and lower groove systems.

The plates’ temperatures are expressed as Fourier expansions of the form

(2.5a)\begin{gather}{\theta _L}(x) = R{a_{uni}} + R{a_{P,L}}{\varTheta _L}(x) = R{a_{uni}} + R{a_{P,L}}\sum\limits_{n ={-} {N_L},\;n \ne 0}^{n ={+} {N_L}} {\theta _L^{(n)}\,{\textrm{e}^{\textrm{i}n(\alpha x + {\varOmega _{T,L}})}}} ,\end{gather}

(2.5b)\begin{gather}{\theta _U}(x) = R{a_{P,U}}{\varTheta _U}(x) = R{a_{P,U}}\sum\limits_{n ={-} {N_L},n \ne 0}^{n ={+} {N_L}} {\theta _U^{(n)}\,{\textrm{e}^{\textrm{i}n(\alpha x + {\varOmega _{T,U}})}}} ,\end{gather}

where ${\varTheta _L}(x)$ and ${\varTheta _U}(x)$ are the shape functions describing temperature distributions at the lower and upper plates, respectively, satisfying conditions

(2.5c)\begin{equation}\textrm{max}[{\varTheta _L}(x)] - \textrm{min}[{\varTheta _L}(x)] = 1,\quad \textrm{max}[{\varTheta _U}(x)] - \textrm{min}[{\varTheta _U}(x)] = 1,\end{equation}

$\theta _L^{(n)}$ and $\theta _U^{(n)}$ are the coefficients of the Fourier expansions describing the temperature profile at the lower and upper plates, respectively, ${N_L}$ is the number of Fourier modes required to describe the form of heating, $R{a_{uni}} = |{g}|\beta {h^3}{\theta _{uni}}/(\kappa \nu )$ is the uniform Rayleigh number measuring the intensity of the uniform component of heating, $R{a_{P,L}} = |{g}|\beta {h^3}{\theta _{LP}}/(\kappa \nu )$ is the lower periodic Rayleigh number measuring the intensity of the lower periodic heating, $R{a_{P,U}} = |{g}|\beta {h^3}{\theta _{UP}}/(\kappa \nu )$ is the upper periodic Rayleigh number measuring the intensity of the upper periodic heating, ${\theta _{LP}}$ and ${\theta _{UP}}$ are the differences between the maximum and minimum of the lower and upper periodic temperature components, ${\varOmega _{T,L}}$ stands for the phase shift between the lower groove and the lower temperature distribution and ${\varOmega _{T,U}}$ stands for the phase shift between the lower groove and the upper-temperature distribution. We focus on systems with perfectly tuned heating and groove modulations; these modulations are described by the same wavenumber $\alpha $.

The stability problem represents an initial-value problem that can be solved using DNS for any initial disturbance velocity and temperature fields set. The Fourier–Chebyshev expansions combined with IBC method provide the required high-accuracy discretization capable of handling complex geometry (Panday & Floryan Reference Panday and Floryan2020, Reference Panday and Floryan2021). The spectral element method described in Cantwell et al. (Reference Cantwell2015) represents an alternative, but its existing implementations provide only a low-order temporal discretization. The DNS requires the specification of a computational box containing an integer number of modulation wavelengths and an integer number of disturbance wavelengths. Such a box may have to be very long. Its size limits possible investigations to a few combinations of groove and heating wavelengths (Blancher et al. Reference Blancher, Le Guer and El Omari2015; Gepner & Floryan Reference Gepner and Floryan2016, Reference Gepner and Floryan2023) and prevents analysis of arbitrary disturbances. Here, we pursue a strategy that avoids these difficulties. We start with determining stationary states and follow up with their stability analysis to determine onset conditions for secondary states.

The following section briefly discusses the algorithm required for determining stationary states. This algorithm has been described elsewhere (Panday & Floryan Reference Panday and Floryan2020, Reference Panday and Floryan2021), so the presentation is limited to the basic concepts required to describe the stability algorithm.

3. Determination of stationary states

Stationary states are described by the stationary version of (2.1) and (2.2). It is helpful to restate these equations for clarity of the presentation:

(3.1a–c)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }{\bar{V}_B} = 0,\quad ({\bar{V}_B}\boldsymbol{\cdot }\boldsymbol{\nabla }){\bar{V}_B} ={-} \boldsymbol{\nabla }{p_B} + {\nabla ^2}{\bar{V}_B} + P{r^{ - 1}}{\theta _B}\bar{j},\quad ({\bar{V}_B}\boldsymbol{\cdot }\boldsymbol{\nabla }){\theta _B} = P{r^{ - 1}}{\nabla ^2}{\theta _B},\end{gather}

(3.2a–d)\begin{gather}{\bar{V}_B}[{y_L}(x)] = {\bar{V}_B}[{y_U}(x)] = 0,\quad {\theta _B}[{y_L}(x)] = {\theta _L}(x),\; \quad {\theta _B}[{y_U}(x)] = {\theta _U}(x),\end{gather}

where subscript B denotes stationary quantities. One needs to add either the pressure gradient constraint ${\wp _x}$ or the flow rate constraint ${\mathrm{\mathbb{Q}}_x}$ for the x direction, and either the pressure constraint ${\wp _z}$ or the flow rate constraint ${\mathrm{\mathbb{Q}}_z}$ for the z direction. Most of the computations have been carried out with ${\wp _x} = 0$ and ${\wp _z} ={-} 2Re$, where the Reynolds number is defined as $Re = {({w_B})_{max}}h/\nu $.

We use spectrally accurate spatial discretization capable of reducing the discretization error to machine accuracy. This discretization relies on Fourier expansions in the z-streamwise and x-spanwise directions and Chebyshev expansions in the y-transverse direction. The algebraic equations are constructed using the Galerkin projection method. The irregularity of the flow domain is handled using the spectrally accurate IBC method, with flow boundary conditions replaced by constraints implemented using the tau method. Details of the algorithm are given in Panday & Floryan (Reference Panday and Floryan2020).

For simplicity of presentation of the stability algorithm, we express velocity components, pressure and temperature of stationary states in the following form:

(3.3)\begin{equation}[{u_B},{v_B},{w_B},{p_B},{\theta _B}](x,y) = \sum\limits_{n ={-} {N_B}}^{{N_B}} {[f_u^{\langle n\rangle },f_v^{\langle n\rangle },f_w^{\langle n\rangle },f_p^{\langle n\rangle },f_\theta ^{\langle n\rangle }](y)\,{\textrm{e}^{\textrm{i}\,n\alpha x}}} .\end{equation}

We have also defined the vorticity vector as

(3.4)\begin{equation}{\bar{\omega }_B} = ({\xi _B},{\eta _B},{\phi _B}) = \left( {\frac{{\partial {w_B}}}{{\partial y}} - \frac{{\partial {v_B}}}{{\partial z}},\frac{{\partial {u_B}}}{{\partial z}} - \frac{{\partial {w_B}}}{{\partial x}},\frac{{\partial {v_B}}}{{\partial x}} - \frac{{\partial {u_B}}}{{\partial y}}} \right),\end{equation}

and write its components as Fourier expansions of the form

(3.5)\begin{equation}[{\xi _{B,\; }}{\eta _B},{\phi _B}](x,y) = \sum\limits_{n ={-} {N_B}}^{{N_B}} {[f_\xi ^{\langle n\rangle },f_\eta ^{\langle n\rangle },f_\phi ^{\langle n\rangle }](y)\,{\textrm{e}^{\textrm{i}\,n\alpha x}}} .\end{equation}

4. Linear stability analysis

The stability analysis begins with the governing equations expressed in terms of the vorticity transport, energy and continuity equations in the following form:

(4.1a,b)\begin{gather}\frac{{\partial \bar{\omega }}}{{\partial t}} - (\bar{\omega }\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{V} + (\bar{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\bar{\omega } = {\nabla ^2}\bar{\omega } + \boldsymbol{\nabla } \times (P{r^{ - 1}}\theta \bar{j}),\quad \frac{{\partial \theta }}{{\partial t}} + (\bar{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\theta = P{r^{ - 1}}{\nabla ^2}\theta ,\end{gather}

(4.1c,d)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\bar{V} = 0,\quad \bar{\omega } = \boldsymbol{\nabla } \times \bar{V}.\end{gather}

Unsteady three-dimensional disturbances are superposed on the stationary flow in the form

(4.2a–c)\begin{align} \bar{\omega } &= {\bar{\omega }_B}(x,y) + {\bar{\omega }_D}(x,y,z,t),\quad \bar{V} = {\bar{V}_B}(x,y) + {\bar{V}_D}(x,y,z,t),\nonumber\\ \theta &= {\theta _B}(x,y) + {\theta _D}(x,y,z,t),\end{align}

where subscript D denotes disturbance quantities. The flow quantities (4.2) are substituted into (4.1), the mean parts are subtracted and the equations are linearized, yielding the linear disturbance equations of the form

(4.3a)\begin{align} &\frac{{\partial {{\bar{\omega }}_D}}}{{\partial t}} - ({\bar{\omega }_B}\boldsymbol{\cdot }\boldsymbol{\nabla }){\bar{V}_D} - ({\bar{\omega }_D}\boldsymbol{\cdot }\boldsymbol{\nabla }){\bar{V}_B} + ({\bar{V}_B}\boldsymbol{\cdot }\boldsymbol{\nabla }){\bar{\omega }_D} + ({\bar{V}_D}\boldsymbol{\cdot }\boldsymbol{\nabla }){\bar{\omega }_B} \nonumber\\ &\quad = {\nabla ^2}{\bar{\omega }_D} + \boldsymbol{\nabla } \times (P{r^{ - 1}}{\theta _D}\bar{j}),\end{align}

(4.3b–d)\begin{align}&\frac{{\partial {\theta _D}}}{{\partial t}} + ({\bar{V}_B}\boldsymbol{\cdot }\boldsymbol{\nabla }){\theta _D} + ({\bar{V}_D}\boldsymbol{\cdot }\boldsymbol{\nabla }){\theta _B} = P{r^{ - 1}}{\nabla ^2}{\theta _D},\quad \boldsymbol{\nabla }\boldsymbol{\cdot }{\bar{V}_D} = 0,\quad {\bar{\omega }_D} = \boldsymbol{\nabla } \times {\bar{V}_D}.\end{align}

The homogeneous boundary conditions of the form

(4.4a,b)\begin{equation}{\bar{V}_D}(x,y,z,t) = 0,\quad {\theta _D}(x,y,z,t) = 0\quad \textrm{at}\;y = {y_L}(x)\textrm{ and }y = {y_U}(x)\end{equation}

complete the formulation.

The solution of (4.3) can be written as a travelling wave with modulated amplitude, i.e.

(4.5a–c)\begin{equation}[{\bar{V}_D},{\bar{\omega }_D},{\theta _D}](x,y,z,t) = [{\bar{G}_D},{\bar{\varOmega }_D},{\kappa _D}](x,y)\,{\textrm{e}^{\textrm{i}(\delta x + \mu z - \sigma t)}} + \textrm{c}\textrm{.c}\textrm{.},\end{equation}

where δ and μ are the real wavenumbers in the x and z directions, respectively, representing components of the disturbance wavevector, $\sigma = {\sigma _r} + \textrm{i}{\sigma _i}$ is the complex frequency with ${\sigma _i}$ describing the rate of growth of disturbances and ${\sigma _r}$ describing their frequency and c.c. stands for the complex conjugate. Functions ${\bar{G}_D}(x,y)$, ${\bar{\varOmega }_D}(x,y)$ and ${\kappa _D}(x,y)$ are the x-periodic amplitude functions accounting for the groove and heating-induced modulations. Substituting (4.5) into (4.3) leads to partial differential equations for the amplitude functions. Since these functions are x-periodic, they can be expressed in terms of the Fourier expansions of the form

(4.6a)\begin{gather}{\bar{G}_D}(x,y) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_u^{\langle m\rangle }(y),g_v^{\langle m\rangle }(y),g_w^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}\,m\alpha x}}} + \textrm{c}\textrm{.c}\textrm{.},\end{gather}

(4.6b)\begin{gather}{\bar{\varOmega }_D}(x,y) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_\xi ^{\langle m\rangle }(y),g_\eta ^{\langle m\rangle }(y),g_\phi ^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}\,m\alpha x}}} + \textrm{c}\textrm{.c}\textrm{.},\end{gather}

(4.6c)\begin{gather}{\kappa _D}(x,y) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_\theta ^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}\,m\alpha x}}} + \textrm{c}\textrm{.c}\textrm{.},\end{gather}

where $g_u^{\langle m\rangle }(y),g_v^{\langle m\rangle }(y),g_w^{\langle m\rangle }(y),g_\xi ^{\langle m\rangle }(y),g_\eta ^{\langle m\rangle }(y),g_\phi ^{\langle m\rangle }(y),g_\theta ^{\langle m\rangle }(y)$ stand for the modal functions.

A combination of (4.5) and (4.6) provides the final expressions for ${{\bar{V}}_D}$, ${{\bar{\omega }}_D}$ and ${{\theta }_D}$ of the form

(4.7a)\begin{gather}{\bar{V}_D}(x,y,z,t) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_u^{\langle m\rangle }(y),g_v^{\langle m\rangle }(y),g_w^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} + \textrm{c}\textrm{.c}\textrm{.},\end{gather}

(4.7b)\begin{gather}{\bar{\omega }_D}(x,y,z,t) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_\xi ^{\langle m\rangle }(y),g_\eta ^{\langle m\rangle }(y),g_\phi ^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} + \textrm{c}\textrm{.c}\textrm{.},\end{gather}

(4.7c)\begin{gather}{\theta _D}(x,y,z,t) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_\theta ^{\langle m\rangle }(y)]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} + \textrm{c}\textrm{.c}\textrm{.}\end{gather}

The above formulation provides means for analysis of the effects of continuous variations of the disturbance wavenumbers δ and μ for a specific modulation wavenumber $\alpha $, i.e. requires a computational box of the same size as the box used to determine the stationary state.

Substitution of (4.7) into (4.3) leads to $(2{N_D} + 1)$ ordinary differential equations (ODEs) for the unknown modal functions for each partial differential equation in the above system. Separation of Fourier modes and extensive rearrangements provide an explicit form of these ODEs expressed in terms of $g_v^{\langle m\rangle },g_\eta ^{\langle m\rangle },g_\theta ^{\langle m\rangle }$ suitable for computations, i.e.

(4.8a)\begin{gather}{T^{\langle m\rangle }}g_v^{\langle m\rangle } - P{r^{ - 1}}k_m^2g_\theta ^{\langle m\rangle } = -\sum\limits_{n ={-} {N_M}}^{{N_M}} [(T_1^{\langle m - n\rangle } + T_2^{\langle m - n\rangle } + T_3^{\langle m - n\rangle })g_v^{\langle m - n\rangle }\nonumber\\ +\,(T_4^{\langle m - n\rangle } - T_5^{\langle m - n\rangle })\; g_\eta ^{\langle m - n\rangle }] ,\end{gather}

(4.8b)\begin{gather}\begin{array}{l} {S^{\langle m\rangle }}g_\eta ^{\langle m\rangle } = \sum\limits_{n ={-} {N_M}}^{ + {N_M}} {[(S_1^{\langle m - n\rangle } + S_2^{\langle m - n\rangle } + S_3^{\langle m - n\rangle })g_v^{\langle m - n\rangle } + (S_4^{\langle m - n\rangle } - S_5^{\langle m - n\rangle })g_\eta ^{\langle m - n\rangle }],} \\ \end{array}\end{gather}

(4.8c)\begin{gather}{Q^{\langle m\rangle }}g_\theta ^{\langle m\rangle } = \sum\limits_{n ={-} {N_M}}^{ + {N_M}} {\left[ {Q_1^{\langle m - n\rangle }g_\theta^{\langle m - n\rangle } + Q_2^{\langle m - n\rangle }g_v^{\langle m - n\rangle } + \frac{{\mu n\alpha f_\theta^{\langle n\rangle }}}{{k_{m - n}^2}}g_\eta^{\langle m - n\rangle }} \right]} ,\end{gather}

where all operators used in the above expressions are defined in Appendix A. The right-hand sides provide modulation-imposed coupling between different modes. The coupling terms contain indices $\langle m - n\rangle$ with the relevant terms truncated for $|m - n|\ge {N_D}$.

The boundary conditions can be written as

(4.9)\begin{equation}\sum\limits_{m ={-} {N_D}}^{{N_D}} {[g_u^{\langle m\rangle },\; g_v^{\langle m\rangle },\; \; g_w^{\langle m\rangle },\; \; g_\theta ^{\langle m\rangle }]\,{\textrm{e}^{\textrm{i}\,m\alpha x}}} = 0\quad \textrm{at}\;y = {y_L}(x)\textrm{ and }y = {y_U}(x).\end{equation}

One needs to extract boundary conditions for different Fourier modes. This process is complex due to the irregularity of the flow domain, which provides coupling between different modes. We provide a detailed explanation in § 4.2.

In the limit of no modulations, the stationary state is reduced to mode zero of ${w_B}$ and ${\theta _B}$. The disturbance wavenumber $\delta $ in the x direction is taken as $\alpha $, and the system is reduced after rearranging indices ${t_m} = (m + 1)\alpha = J\alpha $ to the following form:

(4.10a)\begin{gather}{T^{\langle J\rangle }}g_v^{\langle J\rangle } - P{r^{ - 1}}k_J^2g_\theta ^{\langle J\rangle } + \{ [\textrm{i}\mu {D^2} - \textrm{i}\mu ({D^2} - k_J^2)]f_w^{\left\langle 0 \right\rangle }\} g_v^{\langle J\rangle } = 0,\end{gather}

(4.10b)\begin{gather}{S^{\langle J\rangle }}g_\eta ^{\langle J\rangle } - i\mu f_w^{\langle 0\rangle }g_\eta ^{\langle J\rangle } + \textrm{i}J\alpha Df_w^{\langle 0\rangle }g_v^{\langle J\rangle } = 0,\end{gather}

(4.10c)\begin{gather}{Q^{\langle J\rangle }}g_\theta ^{\langle J\rangle } - \textrm{i}\mu f_w^{\langle 0\rangle }g_\theta ^{\langle J\rangle } - Df_\theta ^{\langle 0\rangle }g_v^{\langle J\rangle } = 0,\end{gather}

where $k_j^2 = {(J\alpha )^2} + {\mu ^2}$ and $|J|= 0,1,2,3, \ldots ,(2{N_D} + 1)$. Equation (4.10a) is the Orr–Sommerfeld operator with thermal coupling term $- P{r^{ - 1}}k_J^2g_\theta ^{\langle J\rangle }$, (4.10b) is the Squire operator, where $\textrm{i}J\alpha Df_w^{\langle 0\rangle }g_v^{\langle J\rangle }$ is the coupling term with the Orr–Sommerfeld operator, and (4.10c) is the energy equation with $Df_\theta ^{\langle 0\rangle }g_v^{\langle J\rangle }$ providing coupling with the Orr–Sommerfeld operator.

4.1. Discretization of the modal equations

Determination of the amplitude functions requires a numerical solution of a system of ODEs described in the previous section. The modal functions are to be represented in terms of Chebyshev expansions. Our preference for the standard definition of the Chebyshev polynomials requires a preliminary mapping of the solution domain $y \in ( - 1 - {y_b},1 + {y_t})$ onto $\hat{y} \in ( - 1,1)$ using transformation $\hat{y} = \varGamma \textrm{[}y - (1 + {y_t})\textrm{]} + 1$ with $\varGamma = 2/(2 + {y_t} + {y_b})$, where ${y_t}$ and ${y_b}$ identify locations of extremities of the upper and lower walls, respectively. The reader may note that this step is not required when boundary modulations are not present. After transformation, the locations of the walls are specified as

(4.11a,b)\begin{align} {\hat{y}_L}(x) &= 1 - \varGamma (2 + {y_t}) + \varGamma {B_L}\sum\limits_{n ={-} {N_G}}^{n = {N_G}} {H_L^{(n)}\,{\textrm{e}^{\textrm{i}n\alpha x}}} ,\nonumber\\ {\hat{y}_U}(x) &= 1 - \varGamma {y_t} + \varGamma {B_U}\sum\limits_{n ={-} {N_G}}^{n = {N_G}} {H_U^{(n)}\,{\textrm{e}^{\textrm{i}n(\alpha x + {\varOmega _G})}}} .\end{align}

It is convenient to cast them in a shorter form for simplicity of the presentation, i.e.

(4.11c)\begin{align}{\hat{y}_L}(x) &= \sum\limits_{n ={-} {N_G}}^{n = {N_G}} {A_L^{(n)}\,{\textrm{e}^{\textrm{i}n\alpha x}}} ,\quad A_L^{(0)} = 1 + \varGamma ( - 2 - {y_t} + {B_L}H_L^{(0)}),\nonumber\\ A_L^{(n)} &= \varGamma {B_L}H_L^{(n)}\quad \textrm{for}\;n \ne 0,\end{align}

(4.11d)\begin{align} {\hat{y}_U}(x) &= \sum\limits_{n ={-} {N_G}}^{n = {N_G}} {A_U^{(n)}\,{\textrm{e}^{\textrm{i}n\alpha x}}} ,\quad \; A_U^{(0)} = 1 + \varGamma ( - {y_t} + {B_U}H_U^{(0)}),\nonumber\\ A_U^{(n)} &= \varGamma {B_U}H_U^{(n)}\quad \textrm{for}\;n \ne 0.\end{align}

Equation (4.8) expressed in terms of $\hat{y}$ assumes the following form:

(4.12a) \begin{gather} {\check{T}{} ^{\langle m\rangle }}g_v^{\langle m\rangle } - P{r^{ - 1}}k_m^2 g_\theta ^{\langle m\rangle } = -\sum\limits_{n ={-}{N_M}}^{{N_M}} [(\check{T}{}_1^{\langle m - n\rangle } + \check{T}{} _2^{\langle m - n\rangle } + \check{T}{} _3^{\langle m - n\rangle }) g_v^{\langle m - n\rangle }\nonumber\\ +\, (\check{T}{} _4^{\langle m - n\rangle } - \check{T}{} _5^{\langle m - n\rangle })\; g_\eta ^{\langle m - n\rangle }] ,\; \end{gather}

(4.12b)\begin{gather}{\check{S}{} ^{\langle m\rangle }}g_\eta ^{\langle m\rangle } = \; \sum \limits_{n ={-} {N_M}}^{ + {N_M}} [(\check{S}{} _1^{\langle m - n\rangle } + \check{S}{} _2^{\langle m - n\rangle } + \check{S}{} _3^{\langle m - n\rangle })g_v^{\langle m - n\rangle } + (\check{S}{} _4^{\langle m - n\rangle } - \check{S}{} _5^{\langle m - n\rangle })g_\eta ^{\langle m - n\rangle }],\end{gather}

(4.12c)\begin{gather}{\check{Q}{} ^{\langle m\rangle }}g_\theta ^{\langle m\rangle } = \sum\limits_{n ={-} {N_M}}^{ + {N_M}} {\left[ {\check{Q}{}_1^{\langle m - n\rangle }g_\theta^{\langle m - n\rangle } + \check{Q}{}_2^{\langle m - n\rangle }g_v^{\langle m - n\rangle } + \frac{{\mu n\alpha f_\theta^{\langle n\rangle }}}{{k_{m - n}^2}}g_\eta^{\langle m - n\rangle }} \right]} ,\end{gather}

where all operators appearing in (4.12) are defined in Appendix B.

System (4.12) has variable coefficients involving stationary state – in the first step of the numerical solution, these coefficients are represented as Chebyshev expansions of the form

(4.13)\begin{equation}[{u_B},{v_B},{w_B}](x,\hat{y}) = \sum\limits_{n ={-} {N_B}}^{{N_B}} {\sum\limits_{r = 0}^{{N_C} - 1} {[G_{u,r}^{\langle n\rangle },\; G_{v,r}^{\langle n\rangle },G_{w,r}^{\langle n\rangle }]{T_r}(\hat{y})\,{\textrm{e}^{\textrm{i}\,n\alpha x}}} } ,\end{equation}

where $G_{u,r}^{\langle n\rangle },\; G_{v,r}^{\langle n\rangle }\textrm{ and }G_{w,r}^{\langle n\rangle }$ are known. In the second step, the unknown modal functions are expressed in terms of the Chebyshev expansions of the form

(4.14)\begin{equation}[g_v^{\langle m\rangle },g_\eta ^{\langle m\rangle },g_\theta ^{\langle m\rangle }](\hat{y}) = \sum\limits_{k = 0}^{{N_T} - 1} {[G_{k,v}^{\langle m\rangle },G_{k,\eta }^{\langle m\rangle },G_{k,\theta }^{\langle m\rangle }]{T_k}(\hat{y})} ,\end{equation}

where $G_{k,v}^{\langle m\rangle },G_{k,\eta }^{\langle m\rangle }\;\textrm{and}\;G_{k,\theta }^{\langle m\rangle }$ are unknown and are to be determined. Substitution of (4.13) and (4.14) into (4.12) and application of the Galerkin projection method result, after a rather lengthy algebra, in a system of algebraic equations for these coefficients of the following form:

(4.15a)\begin{gather}\sum\limits_{k = 0}^{{N_T} - 1} {\{ {A^{\langle m\rangle }}G_{k,v}^{\langle m\rangle } - {H^{\langle m\rangle }}G_{k,\theta }^{\langle m\rangle }\} }\nonumber\\ \quad + \sum\limits_{n ={-} {N_M}}^{{N_M}} {\sum\limits_{k = 0}^{{N_T} - 1} {\sum\limits_{r = 0}^{{N_T} - 1} {\begin{array}{*{20}{c}} {\sum\{ {P^{\langle m - n\rangle }}} \end{array}G_{k,v}^{\langle m - n\rangle } + {R^{\langle m - n\rangle }}G_{k,\eta }^{\langle m - n\rangle }\} } } } = 0,\end{gather}

(4.15b)\begin{gather}\sum\limits_{k = 0}^{{N_T} - 1} {{D^{\langle m\rangle }}G_{k,\eta }^{\langle m\rangle }} + \sum\limits_{n ={-} {N_M}}^{{N_M}} {\sum\limits_{k = 0}^{{N_T} - 1} {\sum\limits_{r = 0}^{{N_T} - 1} {\{ {K^{\langle m - n\rangle }}G_{k,v}^{\langle m - n\rangle } + {Q^{\langle m - n\rangle }}G_{k,\eta }^{\langle m - n\rangle }\} } } } = 0,\end{gather}

(4.15c)\begin{gather} \sum\limits_{k = 0}^{{N_T} - 1} {{I^{\langle m\rangle }}G_{k,\theta }^{\langle m\rangle }} + \sum\limits_{n ={-} {N_M}}^{{N_M}} \sum\limits_{k = 0}^{{N_T} - 1} \sum\limits_{r = 0}^{{N_T} - 1} \{ {U^{\langle m - n\rangle }}G_{k,\theta }^{\langle m - n\rangle } \nonumber\\ +\, {\textrm{V}^{\langle m - n\rangle }}G_{k,v}^{\langle m - n\rangle } + {W^{\langle m - n\rangle }}G_{k,\eta }^{\langle m - n\rangle }\} = 0,\end{gather}

where all operators appearing in (4.15) are defined in Appendix C.

4.2. Discretization of boundary conditions

The discretization of boundary conditions without geometry modulations leads to a standard form of boundary conditions for each modal function. Modulations create coupling between the modal functions – this coupling is different from the coupling provided by the modulations of the field equations discussed previously. The form of these conditions is presented below.

Thermal, no-slip and no-penetration boundary conditions (4.4) at the lower wall are restated as

(4.16a)\begin{gather}{\bar{V}_D}(x,{\hat{y}_L}(x),z,t) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_u^{\langle m\rangle }({{\hat{y}}_L}(x)),g_v^{\langle m\rangle }({{\hat{y}}_L}(x)),g_w^{\langle m\rangle }({{\hat{y}}_L}(x))]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} = 0,\end{gather}

(4.16b)\begin{gather}{\theta _D}(x,{\hat{y}_L}(x),z,t) = \sum\limits_{m ={-} {N_D}}^{ + {N_D}} {[g_\theta ^{\langle m\rangle }({{\hat{y}}_L}(x))]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} = 0.\end{gather}

To be consistent with the modal equations, they must be expressed using the vertical velocity and vorticity components. Their form is given below:

(4.17a)\begin{gather}{u_D}[x,{\hat{y}_L}(x),z,t] = \sum \limits_{m ={-} {N_D}}^{{N_D}} \left[ {\frac{{\textrm{i}\varGamma {t_m}}}{{k_m^2}}g_v^{\langle m\rangle }({{\hat{y}}_L}(x)) - \frac{{\textrm{i}\mu }}{{k_m^2}}g_\eta^{\langle m\rangle }({{\hat{y}}_L}(x))} \right]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}} = 0,\end{gather}

(4.17b)\begin{gather}{v_D}[x,{\hat{y}_L}(x),z,t] = \sum\limits_{m ={-} {N_D}}^{{N_D}} {g_v^{\langle m\rangle }({{\hat{y}}_L}(x))\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} = 0,\end{gather}

(4.17c)\begin{gather}{w_D}[x,{\hat{y}_L}(x),z,t] = \sum \limits_{m ={-} {N_D}}^{{N_D}} \left[ {\frac{{\textrm{i}\mu \varGamma }}{{k_m^2}}g_v^{\langle m\rangle }({{\hat{y}}_L}(x)) + \frac{{\textrm{i}{t_m}}}{{k_m^2}}g_\eta^{\langle m\rangle }({{\hat{y}}_L}(x))} \right]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}} = 0,\end{gather}

(4.17d)\begin{gather}{\theta _D}[x,{\hat{y}_L}(x),z,t] = \sum\limits_{m ={-} {N_D}}^{{N_D}} {g_\theta ^{\langle m\rangle }({{\hat{y}}_L}(x))\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} = 0.\end{gather}

Substitution of the Chebyshev expansion (4.14) into (4.17) results in

(4.18a)\begin{align}&{u_D}[x,{\hat{y}_L}(x),z,t] \nonumber\\ &\quad= \sum\limits_{m ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\varGamma {t_m}}}{{k_m^2}}D{T_k}[{{\hat{y}}_L}(x)]G_{k,v}^{\langle m\rangle } - \frac{{\textrm{i}\mu }}{{k_m^2}}{T_k}[{{\hat{y}}_L}(x)]G_{k,\eta }^{\langle m\rangle }} \right]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} } = 0,\end{align}

(4.18b)\begin{align}&{v_D}[{x,{{\hat{y}}_L}(x ),z,t} ]= \sum \limits_{m ={-} {N_D}}^{{N_D}} \sum \limits_{k = 0}^{{N_T} - 1} {T_k}[{\hat{y}_L}(x )]G_{k,v}^{\left\langle m \right\rangle }{e^{i[{({\delta + m\alpha } )x + \mu z - \sigma t} ]}} = 0, \end{align}

(4.18c)\begin{align}&{w_D}[x,{\hat{y}_L}(x),z,t]\nonumber\\ &\quad = \sum\limits_{m ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\mu \varGamma }}{{k_m^2}}D{T_k}[{{\hat{y}}_L}(x)]G_{k,v}^{\langle m\rangle } + \frac{{\textrm{i}{t_m}}}{{k_m^2}}{T_k}[{{\hat{y}}_L}(x)]G_{k,\eta }^{\langle m\rangle }} \right]\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}}} } = 0,\end{align}

(4.18d)\begin{align}&{\theta _D}[x,{\hat{y}_L}(x),z,t] = \sum \limits_{m ={-} {N_D}}^{{N_D}} \sum \limits_{k = 0}^{{N_T} - 1} {T_k}[{\hat{y}_L}(x)]G_{k,\theta }^{\langle m\rangle }\,{\textrm{e}^{\textrm{i}[(\delta + m\alpha )x + \mu z - \sigma t]}} = 0.\end{align}

The above boundary relations involve values of Chebyshev polynomials and their derivatives at the grooved wall. These values represent periodic functions of the x coordinate and, thus, can be expressed in terms of Fourier expansions as follows:

(4.19a,b)\begin{equation}{T_k}[{\hat{y}_L}(x)] = \sum\limits_{p ={-} {N_j}}^{{N_j}} {({w_L})_k^{\langle p\rangle }\,{\textrm{e}^{\textrm{i}\,p\alpha x}}} ,\quad D{T_k}[{\hat{y}_L}(x)] = \sum\limits_{p ={-} {N_j}}^{{N_j}} {({d_L})_k^{\langle p\rangle }\,{\textrm{e}^{\textrm{i}\,p\alpha x}}} .\end{equation}

Appendix D explains the determination of coefficients ${w_L}$ and ${d_L}$. Substitution of (4.19) into (4.18), extraction of Fourier modes and retention of the first ${N_D}$ modes result in the boundary relations of the form

(4.20a)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\varGamma {t_p}}}{{k_p^2}}({d_L})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle } - \frac{{\textrm{i}\mu }}{{k_p^2}}({w_L})_k^{\langle m - p\rangle }G_{k,\eta }^{\langle p\rangle }} \right]} } = 0,\end{gather}

(4.20b)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {({w_L})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle }} } = 0,\end{gather}

(4.20c)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\mu \varGamma }}{{k_p^2}}({d_L})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle } + \frac{{\textrm{i}{t_p}}}{{k_p^2}}({w_L})_k^{\langle m - p\rangle }G_{k,\eta }^{\langle p\rangle }} \right]} } = 0,\end{gather}

(4.20d)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {({w_L})_k^{\langle m - p\rangle }G_{k,\theta }^{\langle p\rangle }} } = 0,\end{gather}

where truncation $|m - p|\le {N_D}$ is needed to maintain consistency with the modal equations. Increasing the number of retained boundary relations is possible as it improves spatial accuracy but leads to over-constrained formulation (Husain et al. Reference Husain, Szumbarski and Floryan2009). A similar process leads to boundary relations on the upper wall of the form

(4.21a)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\varGamma {t_p}}}{{k_p^2}}({d_U})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle } - \frac{{\textrm{i}\mu }}{{k_p^2}}({w_U})_k^{\langle m - p\rangle }G_{k,\eta }^{\langle p\rangle }} \right]} } = 0,\end{gather}

(4.21b)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {({w_U})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle }} } = 0, \end{gather}

(4.21c)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {\left[ {\frac{{\textrm{i}\mu \varGamma }}{{k_p^2}}({d_U})_k^{\langle m - p\rangle }G_{k,v}^{\langle p\rangle } + \frac{{\textrm{i}{t_p}}}{{k_p^2}}({w_U})_k^{\langle m - p\rangle }G_{k,\eta }^{\langle p\rangle }} \right]} } = 0,\end{gather}

(4.21d)\begin{gather}\sum\limits_{p ={-} {N_D}}^{{N_D}} {\sum\limits_{k = 0}^{{N_T} - 1} {({w_U})_k^{\langle m - p\rangle }G_{k,\theta }^{\langle p\rangle }} } = 0.\end{gather}

Relations (4.20) and (4.21) provide intermodal coupling associated with the grooves.

4.3. The linear algebraic system

We include boundary relations in the system (4.15) using the Tau method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1992), i.e. equations corresponding to the highest-order Chebyshev polynomials in each of (4.15) are replaced with the boundary relations (4.20) and (4.21). The linear system used for computations has a simple form:

(4.22)\begin{equation}\boldsymbol{\varLambda }x = 0,\end{equation}

where $\boldsymbol{\varLambda }$ represents the coefficient matrix and x identifies the vector of unknowns. The structure of the coefficient matrix is illustrated in figure 4(a), where the large blocks identified by green lines correspond to different sets of unknowns (i.e. $v,\eta ,\theta $) – the off-diagonal blocks provide coupling between these unknowns. The structure of a single block for v, displayed in figure 4(b), illustrates couplings between Fourier modes with the diagonal subblocks corresponding to a specific modal function and the off-diagonal subblocks providing couplings between the modes. Rows with green shading identify positions of boundary relations that provide additional intermodal coupling.

Figure 4. (a) Structure of the coefficient matrix $\boldsymbol{\varLambda }$ for ${N_D} = 2$ and ${N_T} = 10$. Green lines identify off-diagonal blocks providing coupling between different unknowns. (b) Structure of a single block with green shading identifying entries corresponding to boundary relations. Black symbols mark the non-zero elements.

5. Method of solution

A homogeneous system (4.22) has a non-trivial solution only for certain combinations of parameters that are interrelated through a dispersion relation of the form

(5.1)\begin{equation}\varTheta \textrm{(}\alpha ,\mu ,\delta ,\sigma ,Re,R{a_{uni}},R{a_{P,L}},R{a_{P,U}},Pr,{B_L},{B_U},H_L^{(n)},H_U^{(n)},\theta _L^{(n)},\theta _U^{(n)}\textrm{)} = 0.\end{equation}

All but two real quantities can be selected arbitrarily. Determining these two quantities requires finding zeros of (5.1). The explicit form of (5.1) is not given. The dispersion relation can be posed in various ways for computational purposes. The first form involves posing the system as an algebraic eigenvalue problem for $\sigma $ in the form

(5.2)\begin{equation}\boldsymbol{AE} = \sigma \boldsymbol{BE},\end{equation}

where $\boldsymbol{E}$ denotes the eigenvector and $\boldsymbol{A}$ and $\boldsymbol{B}$ are the coefficient matrices. The $\sigma $-spectrum can be computed using a general eigenvalue solver (Moler Reference Moler2004). This process can be computationally expensive, and the computed spectrum may suffer from accuracy problems as our stability problem leads to large algebraic systems. The cost can be reduced by evaluating only a spectrum segment using the Arnoldi method (Saad Reference Saad2003). All these methods are suitable for locating eigenvalues of interest but unsuitable for their tracing through parameter space.

Local solutions are more computationally efficient and accurate but typically produce just one eigenvalue (Moler Reference Moler2004). The first method used in this study, the inverse iterations method (Demmel Reference Demmel1997), is suitable for tracing complex frequencies. One starts with an initial approximation for the eigenpair $({\boldsymbol{E}_p},{\sigma _p}) = ({\boldsymbol{E}^{(0)}},{\sigma _0})$ and improve it iteratively. The iteration process uses the following relation:

(5.3)\begin{equation}(\boldsymbol{A} - {\sigma _0}\boldsymbol{B}){\boldsymbol{E}^{(b + 1)}} = \boldsymbol{B}{\boldsymbol{E}^{(b)}},\end{equation}

which starts with $b = 0$ and determines the next approximation for the eigenpair as

(5.4)\begin{equation}\sigma _p^{(b + 1)} = \boldsymbol{E}_p^{(b)T}\boldsymbol{AE}_p^{(b)}/\boldsymbol{E}_p^{(b)\ast }\boldsymbol{BE}_p^{(b)},\end{equation}

where the asterisk denotes the complex conjugate transpose. If ${\sigma _p}$ is the eigenvalue closest to ${\sigma _0}$, then ${\boldsymbol{E}^{(b)}}$ converges to ${\boldsymbol{E}_p}$. The initial approximation ${\boldsymbol{E}^{(0)}}$ must be consistent with boundary conditions.

The Newton–Raphson method posed either in terms of one complex variable or in terms of two real variables offers various alternatives. Search for zeros of the determinant of $\boldsymbol{\varLambda }$ is one of them. This method is ineffective as the determinant varies by several orders of magnitude as a function of the unknown eigenvalue, which leads to numerical difficulties. The alternative approach involves converting the homogeneous problem (4.22) into an inhomogeneous one. To do so, one replaces a homogeneous boundary condition with an inhomogeneous condition, making the system inhomogeneous and easy to solve. The omitted homogeneous boundary condition provides a test verifying if the correct eigenvalue has been selected. This process involves a search for zero of the omitted boundary condition. A simple example of this process in the case of a smooth surface consists of replacing the homogeneous boundary condition for $D{v^{(1)}} = 0$ at one of the walls with an inhomogeneous condition ${D^2}{v^{(1)}} = 1$, which, at the same time, imposes the normalization condition. In general, $D{v^{(1)}}$ is not zero, so one must implement the Newton–Raphson search process to find eigenvalues that bring this condition to zero. Any boundary relation can be used as the test condition. A good initial guess is required to achieve convergence.

6. Testing of the algorithm

The main feature of the algorithm is its high (spectral) accuracy, ability to handle arbitrary temperature and groove patterns very efficiently and ability to handle arbitrary disturbance wavenumbers for any modulation wavenumber. The stationary states were determined using the algorithm described in Panday & Floryan (Reference Panday and Floryan2020). In the tests, the stationary solutions were determined with machine accuracy to eliminate any numerical error in the determinations of these solutions, which may affect the accuracy of stability results.

All the tests reported here have been carried out with zero-pressure-gradient constraints in the x direction. The surface and thermal patterns were of the form

(6.1)\begin{align}\left. {\begin{array}{@{}c@{}} {{y_L} ={-} 1 + \dfrac{{{B_L}}}{2}\textrm{cos}(\alpha x),\quad {y_U} ={+} 1,\quad {\theta_L} = R{a_{uni}} + \dfrac{{R{a_{P,L}}}}{2}\textrm{cos}(\alpha x + {\varOmega_{T,L}})\quad \textrm{and}}\\ {{\theta_U} = \dfrac{{R{a_{P,U}}}}{2}\textrm{cos(}\alpha x + {\varOmega_{T,U}}\textrm{)}\textrm{.}} \end{array}} \right\}\end{align}

The tests dealt with (i) two-dimensional disturbances $(\delta = 0)$ and (ii) three-dimensional disturbances $(\delta \ne 0,\; \mu \ne 0)$.

The algorithm reproduced results for special limiting cases available in the literature (Rayleigh Reference Rayleigh1916; Orszag Reference Orszag1971; Schmid & Henningson Reference Schmid and Henningson2001; Hossain & Floryan Reference Hossain and Floryan2013; Moradi & Floryan Reference Moradi and Floryan2014, Reference Moradi and Floryan2019; Mohammadi et al. Reference Mohammadi, Moradi and Floryan2015).

We provide detailed testing and demonstration of the accuracy of evaluating stability characteristics for flows through a slot formed by grooved walls exposed to periodic heating. The typical spectrum is presented in figure 5, with its width depending on the number of Fourier modes $({N_D})$ used in the computation. There is one unstable eigenvalue for these conditions $({\sigma _r} + \textrm{i}{\sigma _i}) = (458.6385,\; 7.6468)$, which connects to the Squire spectrum in the limit of no modulations.

Figure 5. (a) Spectrum for flow in a channel with longitudinal grooves exposed to spanwise periodic heating with $\alpha = 1$, $\delta = 0.8$, $\mu = 0.5$, $R{a_{P,L}} = 500$, ${\varOmega _{T,L}} = {\rm \pi}/2$, ${B_L} = 0.1$, ${B_U} = 0$, $Re = 1000$. Here ${N_D} = 10$ Fourier modes and ${N_T = 120}$ Chebyshev polynomials were used in the test. The leading mode is identified with a square box. (b) Variations of ${\sigma _i}$ of the leading mode as a function of $R{a_{P,L}}$ for grooved surfaces with ${B_L} = 0.2$ (black colour) with $R{a_{P,L}}$ decreasing to zero and then as a function of ${B_L}$ for isothermal surface (blue colour) with ${B_L}$ decreasing to zero.

The accuracy of the determination of the unstable eigenvalue depends on the number of Fourier modes and Chebyshev polynomials used in the computations. In this testing, the stationary state was determined with machine accuracy, so its accuracy does not affect the stability results. We define error $\mathrm{\Delta }{\sigma _i}$ in the determination of the eigenvalue as

(6.2)\begin{equation}\mathrm{\Delta }{\sigma _i} = |{\sigma _i} - {\sigma _{i,ref}}|,\end{equation}

where ${\sigma _i}$ stands for the computed imaginary part of the eigenvalue and ${\sigma _{i,ref}}$ is its actual value. Since the actual eigenvalue is not known, ${\sigma _{i,ref}}$ has been determined with machine accuracy, which, for the conditions in this test, required the use of ${N_D} = 25$ Fourier modes and ${N_T} = 80$ Chebyshev polynomials.

Figure 6 illustrates variations of $\mathrm{\Delta }{\sigma _i}$ as a function of the number of Fourier modes ${N_D}$ used in the computations. These tests used ${N_T} = 80$ Chebyshev polynomials to reduce this part of discretization error to the machine level. Results for two-dimensional travelling waves show an exponential decrease of the error as the number of Fourier modes increases, demonstrating the spectral accuracy of the algorithm (figure 6a). The absolute value of the error increases with an increase in the heating intensity; however, the character of its exponential decrease is not affected. A similar exponential reduction of error can be observed for the three-dimensional disturbances in figure 6(b). These results indicate that a near machine accuracy can be achieved using ${N_D} = 15$ Fourier modes.

Figure 6. Variations of the error $\mathrm{\Delta }{\sigma _i}$ in the evaluation of the amplification rate ${\sigma _i}$ as a function of the number of Fourier modes ${N_D}$ used in the computations for (a) a two-dimensional travelling wave disturbance $({\delta = 0},\; {\mu = 0.3})$ for $\alpha = 1$, ${B_L} = 0.07,{B_U} = 0$, $Re = 1000$ and (b) a three-dimensional travelling wave disturbance $(\delta = 0.8,\; \mu = 0.5)$ for $\alpha = 1$, ${B_L} = 0.1,{B_U} = 0,{\varOmega _{T,L}} = {\rm \pi}/2,\delta = 0.8,\mu = 0.5$. The heating conditions used in the tests were $R{a_{uni}} = 0$ and $R{a_{P,L}} = R{a_{P,U}} = 50$ and $R{a_{P,L}} = R{a_{P,U}} = 75$. All results were obtained with ${N_T} = 80$ Chebyshev polynomials.

Results displayed in figure 7 illustrate variations of $\mathrm{\Delta }{\sigma _i}$ as a function of the number ${N_T}$ of Chebyshev polynomials used in the computations. These computations used ${N_D = 20}$ Fourier modes to reduce the error associated with this part of discretization to the machine level. The results demonstrate an exponential decrease in the error as the number of Chebyshev polynomials increases. They also illustrate the need to use more Chebyshev polynomials in the case of three-dimensional waves to achieve the same absolute accuracy as in the case of two-dimensional waves. In most cases, 40 Chebyshev polynomials were sufficient to achieve machine accuracy.

Figure 7. Variation of $\mathrm{\Delta }{\sigma _i}$ as a function of the number of Chebyshev polynomials ${N_T}$ for $\alpha = 1$, ${B_L = 0.07},{B_U = 0},{Re = 500}$ for (a) two-dimensional disturbance with $\delta = 0,\mu = 0.4$ and (b) three-dimensional disturbance with $\delta = 0.4,\mu = 0.4$. The heating conditions used in the tests were $R{a_{uni}} = 0$ and $(R{a_{P,L}},R{a_{P,U}}) = (100,\; 0)$ and $(R{a_{P,L}},R{a_{P,U}}) = (200,\; 200)$. All results were obtained with ${N_D} = 20$ Fourier modes.

The IBC method raises the question of how well the homogeneous boundary conditions are enforced at the grooved walls. Since the expected values of the disturbance velocity components and the disturbance temperature are zero, these quantities evaluated at the boundaries represent the absolute errors. They vary in an oscillatory manner along the wall, as illustrated in figure 8. The relative position of the groove and heating patterns determines the location of the maximum error, but the absolute value of the error is never higher than ${10^{ - 10}}$.

Figure 8. Distributions of ${u_D}$ (dash-dotted line), ${v_D}$ (solid line), ${w_D}$ (dashed line) and ${\theta _D}$ (dotted line) at the lower wall for $\alpha = 1$, ${B_L} = 0.07$, ${B_U} = 0$, $Re = 1000$, $R{a_{P,L}} = 100$, $R{a_{P,U}} = 0$, $R{a_{uni}} = 0$, $\mu = 0.4$ and (a) $\delta = 0$ and (b) $\delta = 0.4$. All results were obtained with ${N_D} = 15$ Fourier modes and ${N_T} = 40$ Chebyshev polynomials.

To quantify the boundary error and its variations as a function of the number of Fourier modes, we measured this error using the norm defined as

(6.3)\begin{equation}\parallel\!BE\!\parallel\, = \textrm{max}(q){|_{{y_L}}},\end{equation}

where q is the flow quantity of interest. Eigenfunctions used in (6.3) were normalized by bringing their Euclidean norm defined as

(6.4)\begin{equation}\Vert\boldsymbol{E}\Vert_{2} = \sqrt {\boldsymbol{E}_1^2 + \boldsymbol{E}_2^2 + \cdots + \boldsymbol{E}_{3{N_T}(2{N_D} + 1)}^2} = 1\end{equation}

to unity. Here $\boldsymbol{E}$ is an eigenvector that includes each disturbance quantity appearing in (6.3) with length $3{N_T}(2{N_D} + 1)$. Results displayed in figure 9 demonstrate that ${\parallel} BE\parallel $ decreases exponentially with ${N_D}$. The convergence rate varies between different disturbance quantities, but all have exponential forms.

Figure 9. Variations of the boundary errors as functions of the number ${N_D}$ of Fourier modes used in the computations for $\alpha = 1$, ${B_L} = 0.07$, $Re = 1000$, $R{a_{P,L}} = 100$, $\mu = 0.4$ and (a) $\delta = 0$ and (b) $\delta = 0.4$. All results were obtained with ${N_T} = 80$ Chebyshev polynomials.

Another way to demonstrate the spectral convergence of the algorithm is to determine variations of the Chebyshev norm defined as

(6.5)\begin{equation}\Vert \varPhi _q^{(m)}(\hat{y})\Vert = {\left\{ {\int_{ - 1}^1 {G_{k,q}^{(m)}(\hat{y})\; \overline {G_{k,q}^{( - m)}(\hat{y})} \omega (\hat{y})\,\textrm{d}\hat{y}} } \right\}^{1/2}},\end{equation}

where q stands for the modal function of choice. Figure 10 illustrates variations of this norm with an increase of the mode number m for two-dimensional and three-dimensional disturbances. These results confirm the exponential reduction of the Chebyshev norm as the mode number increases for all the quantities.

Figure 10. Variations of the Chebyshev norms as functions of the mode number for $\alpha = 1$, ${B_L} = 0.07$, ${B_U} = 0$, $Re = 1000$, $R{a_{P,L}} = 100$, $R{a_{P,U}} = R{a_{uni}} = 0$, $\mu = 0.4$ and (a) $\delta = 0$ and (b) $\delta = 0.4$. All results were obtained with ${N_D} = 15$ Fourier modes and ${N_T} = 80$ Chebyshev polynomials.

All the tests reported in this section validate the spectral accuracy of the algorithm. The algorithm has a significant advantage in efficiently handling variations of geometry and heating patterns compared with the grid-based approaches. Its gridless nature makes it suitable for a parametric study involving the analysis of any parameter that defines geometric and heating patterns. Mesh construction and grid independence studies required by analyses using the conventional grid methods are very labour- and time-expensive.

Tables 1 and 2 provide numerical values illustrating the convergence of the test eigenvalue determined using different numbers of Chebyshev polynomials and Fourier modes. The comparison dataset uses a form similar to that of Orszag (Reference Orszag1971), which others can use for verification of the accuracy of their results.

Table 1. Variations of the leading eigenvalue $\sigma $ as a function of the number of Fourier modes ${N_M}$ and the number of Chebyshev polynomials ${N_T}$ used in the computations for an isothermal system with $R{a_{pL}} = 0$, ${B_L} = 0.1$, ${B_U} = 0$, $\alpha = 1$, $\delta = 0$ and $Re = 3000$, $\mu = 0.4$, and $Re = 7000$, $\mu = 1$.

Table 2. Variations of the leading eigenvalue $\sigma $ as a function of the number of Fourier modes ${N_M}$ and the number of Chebyshev polynomials ${N_T}$ used in the computations for a non-isothermal system with $Re = 1000$, $R{a_{pL}} = 500$, ${B_L} = 0.1$, ${B_U} = 0$, $\alpha = 1$, $\delta = 0.8$.

7. Example of an application problem

This section discusses an application problem selected to illustrate the capabilities of the proposed algorithm. This problem cannot be convincingly analysed using the DNS-type methodology.

Consider isothermal pressure-gradient-driven flow in a channel modified by longitudinal grooves. This flow is subject to a new type of instability, which, for simplicity, we refer to as Floryan's instability (Mohammadi et al. Reference Mohammadi, Moradi and Floryan2015). The presence of longitudinal grooves creates spanwise gradients of the streamwise velocity, which activate inflection-point-type inviscid instability. The critical Reynolds number can be as low as O(10) if grooves with large enough amplitudes are used (Gepner & Floryan Reference Gepner and Floryan2020). The instability creates waves with phase speed near the maximum flow velocity, leading to laminar chaos in the saturation state (Gepner & Floryan Reference Gepner and Floryan2020). Use of grooves with amplitude ${B_L} = 0.1$ results in a critical Reynolds number of about $2600$ (see figure 11a). If one adds heating with amplitude $R{a_{P,L}} = 700$ and wavenumber matching the groove wavenumber, and places this heating in such a manner so that hot spots overlap with the groove peaks, e.g. ${\varOmega _{T,L}} = 0$, the critical Reynolds number is reduced to about 375 (see figure 11b). It is thus possible to achieve low critical Re with much smaller groove amplitudes. The use of heating patterns is attractive as heating can be turned on and off as required, and its implementation is likely more economically effective as the creation of special conduit geometries is avoided. Figure 12(a) displays the spectrum for $Re = 400$, which shows the presence of a single unstable eigenvalue. Figure 12(b) illustrates variations in the amplification rate ${\sigma _i}$ as a function of the streamwise wavenumber $\delta $, which is essential for determining $\delta $ giving the largest amplification, this step being required in a stability analysis. Table 3 gives spanwise size W of the computational domain that needs to be used for DNS to reproduce some of these results – this table illustrates difficulties associated with using DNS.

Figure 11. The neutral curves in the $(Re,\mu )$ plane for flow in (a) an isothermal $(R{a_{P,L}} = R{a_{P,U}} = 0)$ grooved ($\alpha = 1$, ${B_L} = 0.1$, ${B_U} = 0$) channel for $\delta = 1$ and (b) the same channel heated from below ($R{a_{P,L}} = 700$, $R{a_{P,U}} = 0$).

Figure 12. (a) Spectrum for flow in a grooved channel exposed to periodic heating for $\alpha = 1$, $R{a_{P,L}} = 700$, $Re = 400$, ${\varOmega _{T,L}} = 0$, ${B_L} = 0.1$, ${B_U} = 0$ and $\delta = 1$, $\mu = 1$ determined using ${N_D} = 10$ Fourier modes and ${N_T} = 120$ Chebyshev polynomials. (b) Variation of the amplification rate ${\sigma _i}$ as a function of $\delta $ for $\mu = 1$ with other conditions as in (a).

Table 3. Spanwise extent W of the computational domain required for the use of the DNS approach for the heating wavenumber $\alpha = 1$ and selected disturbance wavenumbers $\delta $.

8. Summary

An algorithm suitable for analysing the stability of spatially modulated shear layers has been proposed and described in the context of modulations induced by surface topography and modulations created by heating patterns. The former represents geometry modulations, while the latter represents field modulations. The algorithm can be easily extended to other forms of modulations. Spatial modulation problems are characterized by the formation of commensurate and incommensurate states arising out of the need to analyse all possible disturbances, and these states are not accessible to classical DNS-based approaches. The proposed algorithm bypasses this problem by treating disturbances as modulated waves and focusing on determining the amplitude functions. This determination relies on spatial discretization based on Fourier expansions in the streamwise and spanwise directions and Chebyshev expansions in the wall-normal direction. The Galerkin projection method has been used to develop the linear system of algebraic equations for the expansion coefficients. The homogeneous boundary conditions to be imposed along the corrugated walls have been replaced by constraints (the IBC method) included in the complete system of linear equations using the Tau concept. The resulting eigenvalue problem was solved using standard methods. Numerous tests demonstrated that the algorithm provides spectral accuracy. Detailed comparison tables for selected test cases have been provided, and the eigenvalue accuracy has been demonstrated. An example of a physical problem was presented to illustrate the power and effectiveness of the proposed method.

The proposed algorithm is gridless and requires minimal user time to adapt to new geometries of the bounding walls and new heating patterns. It bypasses the need for cumbersome grid convergence studies to verify grid-based methods’ accuracy. The number of Fourier modes and Chebyshev polynomials used in the computations sets the absolute accuracy.