1 Introduction

The pioneering experiments of Schubauer and Skramstad [1] demonstrated for the first time that two-dimensional boundary layers were unstable to a type of instability now known as Tollmien–Schlichting instability. Following their work, there have been numerous experimental, theoretical and computational studies looking at the subsequent non-linear development of the instability and the role it plays in the transition to turbulence, see for instance [2,3,4]. The success of the experimental work [1] in being able to observe instability waves triggered by disturbances generated by a vibrating ribbon was due in part to using a low turbulence wind tunnel. This meant that instability waves which were amplified downstream were picked up by the instruments and were seen to be directly responsible for causing transition to turbulence. In fully three-dimensional boundary layer flows the route to transition can be different as discussed in [5, 6].

Disturbances can be triggered through many different mechanisms such as surface heating [7], free-stream turbulence [8], wall roughness [9], leading edges [10], elastic vibrations [11] and not just a vibrating ribbon. The experiments of [7, 12] showed also how surface heating can be used not only to generate and amplify Tollmien–Schlichting waves, but also control the instability. In these experiments two heating strips were placed on a flat surface and Tollmien–Schlichting waves excited at the first strip were either amplified or substantially reduced in amplitude by an active control applied to the second heating strip. Flow control seeks not only to control instabilities but also change flow properties such as separation and drag, which can in turn lead to more favourable flow for many applications, such as the design of laminar flow wings. The extensive survey of Löfdahl and Gad-el-Hak [13] has reviewed many of the technological and other practical challenges in this area.

One of the objectives of this paper is to study how localised surface heating can be used to control the growth of three-dimensional Tollmien–Schlichting disturbances. The recent paper [14] investigated how surface heating can be used to control and cancel out two-dimensional Tollmien–Schlichting waves excited in a flat plate boundary layer flow. The approach was based on asymptotic methods combined with triple-deck scalings to model a two-dimensional hump-shaped vibrator placed on a flat plate, together with a localised surface heating strip. This was an extension of the problem first looked at by Terent’ev [15, 16] but with thermal effects present. The mathematical analysis presented in [15, 16] showed how unstable spatially growing Tollmien–Schlichting waves could be generated by the vibrator and the author was able to use for his analysis the linearised stability results of Smith [17] and the classical results due to Lin [18], amongst others. The results of [14] showed that it is possible to choose heating profiles such that no Tollmien–Schlichting wave is generated, and an exact expression for the appropriate choice of the heating function stemming from the linear analysis was derived. They proposed an approximate form of the heating profile which was used successfully in numerical computations to significantly reduce unstable Tollmien–Schlichting wave amplitudes. The papers of Terent’ev [15, 16], Brennan et al. [14] are for two-dimensional disturbances in a two-dimensional boundary layer flow. In the current work we study an otherwise two-dimensional compressible boundary layer flow encountering a three-dimensional vibrator and with localised surface heating present. We adopt the same methodology and scalings as in [14, 16] which then leads to the three-dimensional unsteady triple-deck equations governing the linear and non-linear development of the three-dimensional disturbances. The corresponding three-dimensional analogue, including the steady-state version, of the Terent’ev [15, 16] problem has been studied before by [19,20,21,22], and others including [23] but without thermal effects. In the current paper thermal effects are included and we derive a three-dimensional extension of the cancellation formula obtained in Brennan et al. [14] before presenting numerical solutions of the linearised unsteady triple-deck equations using an approximate unsteady heating distribution based on the exact cancellation function. This is shown to work effectively in substantially reducing unstable three-dimensional Tollmien–Schlichting wave amplitudes for several cases.

In Sect. 2 we derive the governing asymptotic unsteady equations valid for large values of the Reynolds number, and investigate the solution properties of the linearised initial-value problem. In Sect. 3 the initial-value problem is solved numerically and results are presented for a range of vibration frequencies and subsonic Mach numbers. Finally in Sect. 4 we finish with brief conclusions.

2 Problem formulation

Consider the laminar subsonic flow past a flat plate containing a vibrator at a distance L from the leading edge of the plate, see Fig. 1. The vibrator is represented by a small patch on the surface of the plate performing oscillations in the normal direction. The plate also contains a localised heating element whose dimensions in the wall normal direction are small compared with the thickness of the boundary layer. We will assume that the Reynolds number Re is large (where \( Re={\rho _{\infty }U_{\infty }L/\mu _{\infty }}\)). At large distances from the plate the flow is uniform with speed \(U_{\infty }\) parallel to the plate, density \(\rho _{\infty }\) and dynamic viscosity coefficient \(\mu _{\infty }\). Upstream of the vibrator the oncoming flow is a two-dimensional boundary layer. The presence of the vibrator and localised heating induces a local three-dimensional interaction and variation of the fluid-dynamic functions in the spanwise direction \(z^*\).

Fig. 1
figure 1

a Schematic of the flow showing the hump-shaped three-dimensional vibrator and the oncoming two-dimensional boundary layer flow \(U_B\) with the triple-deck region. b A side view of the same configuration showing the triple-deck and heating

We non-dimensionalise the variables and flow quantities with respect to a lengthscale L, velocity \(U_{\infty }\), and free-stream density \(\rho _{\infty }\) so that

$$\begin{aligned} x= & {} \frac{x^{*}}{L}, \quad y=\frac{y^{*}}{L},\quad z=\frac{z^{*}}{L}, \quad t=\frac{U_{\infty }}{L}t^*, \\ u= & {} \frac{u^{*}}{U_{\infty }},\quad v=\frac{v^{*}}{U_{\infty }},\quad w=\frac{w^*}{U_{\infty }}, \quad T=\frac{{ T}^{*}{{{\mathcal {R}}}} M_{\infty }^2 \gamma }{U_{\infty }^2},\quad \\ p= & {} \frac{p^{*}-p_{\infty }}{\rho _{_{\infty }}U_{\infty }^{2}},\quad \mu =\frac{\mu ^{*}}{\mu _{\infty }} \quad \mathrm{and} \quad \rho =\frac{\rho ^{*}}{\rho _{\infty }}. \end{aligned}$$

The superscript asterisk quantities are dimensional, (xyz) are the coordinates in the streamwise and wall normal and spanwise directions with corresponding velocity components (uvw),  t is time, T is the temperature, p the pressure, \(\mu \) the dynamic viscosity, \(\rho \) is the density, \(p_{\infty }\) is the free-stream pressure, \({{{\mathcal {R}}}}\) is the specific gas constant, and \(\gamma \) is the ratio of specific heats.

The continuity, Navier–Stokes, energy equation and equation of state are given by

$$\begin{aligned}&{\partial \rho \over \partial t}+{\partial (\rho u)\over \partial x}+{\partial (\rho v)\over \partial y}+ {\partial (\rho w)\over \partial z}=0, \end{aligned}$$
(1a)
$$\begin{aligned}&\rho \left( {\partial u \over \partial t}+u{\partial u\over \partial x}+v{\partial u\over \partial y} +w{\partial u\over \partial z}\right) = -{\partial p\over \partial x} + {1\over Re}{\partial \over \partial y}\left( \mu \left[ {\partial u\over \partial y} +{\partial v\over \partial x}\right] \right) +\cdots \end{aligned}$$
(1b)
$$\begin{aligned}&\rho \left( {\partial v \over \partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right) = -\frac{\partial p}{\partial y}+\cdots \end{aligned}$$
(1c)
$$\begin{aligned}&\rho \left( {\partial w \over \partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right) = -\frac{\partial p}{\partial z}+\frac{1}{Re}\frac{\partial }{\partial y}\left( \mu \left[ \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right] \right) +\cdots , \end{aligned}$$
(1d)
$$\begin{aligned}&\rho \left( {\partial T \over \partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}\right) = ( {\gamma -1})M_{\infty }^2\left( \frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}+v\frac{\partial p}{\partial y}+w\frac{\partial p}{\partial z}\right) +\frac{1}{Pr Re} \frac{\partial }{\partial y}\left( \mu \frac{\partial T}{\partial y}\right) +\cdots \end{aligned}$$
(1e)
$$\begin{aligned}&\gamma M_{\infty }^2 p=\rho T-1. \end{aligned}$$
(1f)

In the above we have only kept the important terms which are relevant to the high-Reynolds number analysis given below and the ‘\(+\cdots \)’ signify terms that can be neglected at leading order.

In addition to the Reynolds number Re we have the Prandtl number Pr,  and \(M_{\infty }\) the free-stream Mach number. Here, \(M_{\infty }= {U_{\infty }}/{c_{\infty }}\) where \(c_{\infty }=\sqrt{\gamma p_{\infty }/\rho _{\infty }}\) is the speed of sound in the undisturbed flow.

As in [14, 24, 25], we take the dimensions of the vibrator and localised heating to be commensurate with the triple-deck scales, such that the vibrator occupies a region of extent \(O(Re^{-3/8})\) in the streamwise and spanwise directions and \(O(Re^{-5/8})\) in the wall normal direction. The reasons for this are explained clearly in the work of Lipatov [25]. The vibrator and localised surface heating induce non-linear interactions and changes to flow properties which can be studied via the triple-deck structure. Other distinguished scales and regimes emerge as asymptotic limits.

We choose the time scale of \(O(Re^{-1/4})\) such that the unsteady terms balance with the convective terms in (1b). Given the earlier comments on the time scale we let \(t=Re^{-{1\over 4}}t_*\) and \(t_*\) is O(1). The unstable and most amplified Tollmien–Schlichting disturbances have frequencies of this order and following the work of [14] this is the most natural starting point.

With the above scalings, the vibrator is situated near \(x=1\) and hence we may set \(x=1+Re^{-{3\over 8}} x_{*}, \quad z=Re^{-{3\over 8}}z_* \) where \(x_*, z_*\) are O(1) in the region occupied by the vibrator.

We will write the equation for the surface of the vibrator as

$$\begin{aligned} y=y_w(t,x,z)=Re^{-{5\over 8}}f(t_*,x_*,z_*) . \end{aligned}$$

We will assume that \(f(t_*,x_*,z_*)\) is zero apart from a finite interval in \(x_*\) and \(z_*\) in the vicinity of \(x=1\), \(z=0\).

The presence of the vibrator means that at the wall we require

$$\begin{aligned} u=w=0, \quad v={\partial y_w \over \partial t} \quad { \mathrm on}\quad y=y_w(t,x,z). \end{aligned}$$

Next we assume that there is also a localised heating element on the flat plate co-located with the vibrator. This is modelled by the wall temperature profile being given by \(T= T_w(t_*,x_*,z_*)\) on \(y=y_w\), where \(T_w\) is an order one quantity. Under these conditions we can use the triple-deck model where the flow is divided into three layers: the near-wall viscous layer (lower deck), the main part of the boundary layer (main deck) and the inviscid region outside the boundary layer (upper deck).

2.1 Main deck

Firstly in the main part of the boundary layer with \(y=Re^{-1/2}y_2\) we have

$$\begin{aligned} u= & {} U_{B}(1,y_{2})+Re^{-{1\over 8}}u_{21}(t_*,x_*,y_2,z_*)+\cdots ,\\ v= & {} Re^{-{1\over 4}}v_{21}(t_*,x_*,y_2,z_*)+\cdots , \\ w= & {} Re^{-{1\over 4}}w_{21}(t_*,x_*,y_2,z_*)+\cdots ,\\ p= & {} Re^{-{1\over 4}}p_{21}(t_*,x_*,y_2,z_*)+\cdots ,\\ \rho= & {} \rho _{B}(1,y_{2})+Re^{-{1\over 8}}\rho _{21}(t_*,x_*,y_2,z_*)+\cdots ,\\ T= & {} T_{B}(1,y_{2})+Re^{-{1\over 8}}T_{21}(t_*,x_*,y_2,z_*)+\cdots . \end{aligned}$$

where the suffix B denotes the oncoming two-dimensional Blasius boundary layer flow. Substitution into the Navier–Stokes equations and solving yields the solutions:

$$\begin{aligned} u_{21}= & {} A_*(t_*,x_*,z_*){\mathrm{d} U_B \over \mathrm{d}y_2}, \quad v_{21}= -{\partial A_* \over \partial x_*} U_B, \quad w_{21}={D_*(t_*,x_*,z_*)\over U_B \rho _B}\nonumber \\ T_{21}= & {} A_* {\mathrm{d}T_B \over \mathrm{d}y_2}, \quad \rho _{21}=A_* {\mathrm{d} \rho _B \over \mathrm{d}y_2},\quad p_{21}=P_*(t_*,x_*,z_*) . \end{aligned}$$
(2)

Here \(A_*(t_*,x_*,z_*), D_*(t_*,x_*,z_*)\) and \(P_*(t_*,x_*,z_*)\) are unknown functions and

$$\begin{aligned} {\partial P_* \over \partial z_*}=-{\partial D_*\over \partial x_*}. \end{aligned}$$

2.2 Upper deck

In the upper deck \(y=Re^{-{3\over 8}}y_{1}\) and the expansions of the fluid-dynamic functions take the form

$$\begin{aligned} u= & {} 1+Re^{-{1\over 4}}u_{11}(t_*,x_*,y_1,z_*)+\cdots ,\\ v= & {} Re^{-{1\over 4}}v_{11}(t_*,x_*,y_1,z_*)+\cdots ,\\ w= & {} Re^{-{1\over 4}}w_{11}(t_*,x_*,y_1,z_*)+\cdots ,\\ p= & {} Re^{-{1\over 4}}p_{11}(t_*,x_*,y_1,z_*)+\cdots ,\\ \rho= & {} 1+Re^{-{1\over 4}}\rho _{11}(t_*,x_*,y_2,z_*)+\cdots ,\\ T= & {} 1+Re^{-{1\over 4}}T_{11}(t_*,x_*,y_1,z_*)+\cdots . \end{aligned}$$

After substituting into the Navier–Stokes equations it is found that the pressure \(p_{11}\) satisfies the equation

$$\begin{aligned} (1-M_{{\infty }}^{2})\frac{\partial ^{2}p_{11}}{\partial x_*^{2}}+\frac{\partial ^{2}p_{11}}{\partial y_{1}^{2}}+\frac{\partial ^{2}p_{11}}{\partial z_*^{2}}=0 \end{aligned}$$

with

$$\begin{aligned}&p_{11}(t_*,x_*,y_1=0,z_*)=P_*(t_*,x_*,z_*),\\&{\partial p_{11}\over \partial y_1}(t_*,x_*,y_1=0,z_*)={\partial ^2 A_* \over \partial x_*^2}, \quad {\partial p_{11}\over \partial z_*}(t_*,x_*,y_1=0,z_*)=-{\partial D_* \over \partial x_*} \end{aligned}$$

arising from matching with the main deck. We also require \(p_{11}(t_*,x_*,y_1,z_*)\rightarrow 0\) as \(x_*^2+y_1^2+z_*^2 \rightarrow \infty \).

2.3 Lower deck

In the lower deck with \(y=Re^{-{5\over 8}}y_3\) the expansions are

$$\begin{aligned} u= & {} Re^{-{1\over 8}}u_{31}(t_*,x_*,y_3,z_*)+\cdots ,\\ v= & {} Re^{-{3\over 8}}v_{31}(t_*,x_*,y_3,z_*)+\cdots ,\\ w= & {} Re^{-{1\over 8}}w_{31}(t_*,x_*,y_3,z_*)+\cdots ,\\ p= & {} Re^{-{1\over 4}}p_{31}(t_*,x_*,y_3,z_*)+\cdots ,\\ \rho= & {} \rho _{31}(t_*,x_*,y_3,z_*)+\cdots ,\\ T= & {} T_{31}(t_*,x_*,y_3,z_*)+\cdots . \end{aligned}$$

Substituting these expansions into (1) we obtain

$$\begin{aligned}&\frac{\partial \rho _{31}}{\partial t_*}+\frac{\partial (u_{31}\rho _{31})}{\partial x_*}+\frac{\partial (v_{31}\rho _{31})}{\partial y_{3}}+\frac{\partial (w_{31}\rho _{31})}{\partial z_*}=0, \end{aligned}$$
(3a)
$$\begin{aligned}&\rho _{31}\left( \frac{\partial u_{31}}{\partial t_*}+u_{31}\frac{\partial u_{31}}{\partial x_*}+v_{31}\frac{\partial u_{31}}{\partial y_{3}}+w_{31}\frac{\partial u_{31}}{\partial z_*}\right) =-\frac{\partial p_{31}}{\partial x_*}+\frac{\partial }{\partial y_{3}}\left( \mu \frac{\partial u_{31}}{\partial y_{3}}\right) , \end{aligned}$$
(3b)
$$\begin{aligned}&\rho _{31}\left( \frac{\partial w_{31}}{\partial t_*}+u_{31}\frac{\partial w_{31}}{\partial x_*}+v_{31}\frac{\partial w_{31}}{\partial y_{3}}+w_{31}\frac{\partial w_{31}}{\partial z_*}\right) =-\frac{\partial p_{31}}{\partial z_*}+\frac{\partial }{\partial y_{3}}\left( \mu \frac{\partial w_{31}}{\partial y_{3}}\right) , \end{aligned}$$
(3c)
$$\begin{aligned}&\frac{\partial p_{31}}{\partial y_{3}}=0, \end{aligned}$$
(3d)
$$\begin{aligned}&\rho _{31}\left( \frac{\partial T_{31}}{\partial t_*}+u_{31}\frac{\partial T_{31}}{\partial x_*}+v_{31}\frac{\partial T_{31}}{\partial y_{3}}+w_{31}\frac{\partial T_{31}}{\partial z_*}\right) =\frac{\partial }{\partial y_{3}}\left( \frac{\mu }{Pr}\frac{\partial T_{31}}{\partial y_{3}}\right) , \end{aligned}$$
(3e)
$$\begin{aligned}&\rho _{31}T_{31}=1 . \end{aligned}$$
(3f)

To match with the main-deck solutions we require that as \(y_3 \rightarrow \infty \),

$$\begin{aligned} \begin{array}{ccc} u_{31}\rightarrow \lambda (A_*(t_*,x_*,z_*)+y_{3}),&T_{31}\rightarrow T_{B0}&\mathrm {and}\,w_{31}\rightarrow \frac{D_*(t_*,x_*,z_*)}{\lambda \rho _{B0}y_{3}},\end{array} \end{aligned}$$

where \(T_{B0}=T_B(1,0), \rho _{B0}=\rho _B(1,0).\) Also as \(x_*\rightarrow -\infty \)

$$\begin{aligned} u_{31} \rightarrow \lambda y_3, \quad A_*(t_*,x_*,z_*) \rightarrow 0,\quad D_*(t_*,x_*,z_*) \rightarrow 0 . \end{aligned}$$

Here \(\lambda = {\partial U_B \over \partial y_3}(y_3=0).\) The additional boundary conditions on the hump-shaped wall are

$$\begin{aligned} u_{31}=w_{31}=0, \quad v_{31}=\frac{\partial f}{\partial t_*}, \quad \mathrm{and}\quad T_{31}=T_w(t_*,x_*,z_*) \quad \mathrm{on} \quad y_3= f(t_*,x_*,z_*) \end{aligned}$$

with \(T_w\) being some prescribed unsteady wall temperature.

We first make use of the combined Prandtl Dorodnitsyn–Howarth transform given by

$$\begin{aligned} \rho _{31}v_{31} = v_* -{\partial y_*\over \partial t_*}-u_{31}{\partial y_* \over \partial x_*} -w_{31}{\partial y_* \over \partial z_*}, \end{aligned}$$

where

$$\begin{aligned} y_* = \int _{f}^{y_3}\rho _{31}(t_*,x_*,y_3,z_*)\,\mathrm{d}y_3 . \end{aligned}$$

The form of the transform above with the lower limit non-zero is a novel extension to the usual Dorodnitsyn–Howarth transform and as used also by [14] where further details may be found. The equations (3) and boundary conditions reduce to

$$\begin{aligned}&{\partial {u}_{31}\over \partial x_*}+ {\partial { v}_*\over \partial {y_*}}+{\partial w_{31} \over \partial z_*}=0, \end{aligned}$$
(4a)
$$\begin{aligned}&{\partial {u}_{31} \over \partial t_*}+ {u}_{31}{\partial {u}_{31}\over \partial x_*}+ {v}_* {\partial {u}_{31}\over \partial {y_*}}+ w_{31} {\partial u_{31} \over \partial z_*}= -T_{31} {\partial {p}_{31}\over \partial x_*} + {\partial \over \partial y_*}\left( \rho _{31}\mu {\partial {u}_{31}\over \partial y_* }\right) , \end{aligned}$$
(4b)
$$\begin{aligned}&0 = -{\partial {p}_{31}\over \partial {y_*}}, \end{aligned}$$
(4c)
$$\begin{aligned}&{\partial {w}_{31} \over \partial t_*}+ {u}_{31}{\partial {w}_{31}\over \partial x_*}+ {v}_* {\partial {w}_{31}\over \partial {y_*}}+ w_{31} {\partial w_{31} \over \partial z_*}= -T_{31}{\partial {p}_{31}\over \partial z_*} + {\partial \over \partial y_*}\left( \rho _{31}\mu {\partial {w}_{31}\over \partial y_* }\right) , \end{aligned}$$
(4d)
$$\begin{aligned}&{\partial {T}_{31} \over \partial t_*}+ {u}_{31}{\partial {T}_{31}\over \partial x_*}+ {v}_* {\partial {T}_{31}\over \partial {y_*}}+ w_{31}{\partial T_{31} \over \partial z_*}={1\over Pr}{\partial \over \partial y_*}\left( \rho _{31}\mu {\partial {T}_{31} \over \partial y_*}\right) , \end{aligned}$$
(4e)
$$\begin{aligned}&\rho _{31}T_{31} = 1 . \end{aligned}$$
(4f)

The boundary conditions are

$$\begin{aligned}&v_*(t_*,x_*,y_*,z_*)=0 \quad \mathrm{on} \quad y_*=0, \end{aligned}$$
(5a)
$$\begin{aligned}&u_{31}(t_*,x_*,y_*,z_*)=0 \quad \mathrm{on }\quad y_*=0, \end{aligned}$$
(5b)
$$\begin{aligned}&w_{31}(t_*,x_*,y_*,z_*)=0 \quad \mathrm{on }\quad y_*=0, \end{aligned}$$
(5c)
$$\begin{aligned}&T_{31}(t_*,x_*,y_*,z_*)=T_w(t_*,x_*,z_*) \quad \mathrm{on }\quad y_*=0, \end{aligned}$$
(5d)
$$\begin{aligned}&T_{31}(t_*,x_*,y_*,z_*)=T_{B0} \quad \mathrm{as }\quad y_*\rightarrow \infty , \end{aligned}$$
(5e)
$$\begin{aligned}&u_{31}(t_*,x_*,y_*,z_*)= \lambda T_{B0} y_* \quad \mathrm{as}\quad x_* \rightarrow -\infty , \end{aligned}$$
(5f)
$$\begin{aligned}&w_{31}(t_*,x_*,y_*,z_*)= {D_* \over \lambda y_*} \quad \mathrm{as}\quad y_* \rightarrow \infty , \end{aligned}$$
(5g)
$$\begin{aligned}&u_{31}(t_*,x_*,y_*,z_*)=\lambda (T_{B0}y_* + K_*(t_*,x_*,z_*)) \quad \mathrm{as}\quad y_* \rightarrow \infty , \end{aligned}$$
(5h)

where

$$\begin{aligned} K_*(t_*,x_*,z_*)= f(t_*,x_*,z_*)+\int _0^{\infty }( T_{31} -T_{B0})\,\mathrm{d}y_* + A_*(t_*,x_*,z_*) . \end{aligned}$$
(5i)

The term involving the integral in the expression for \(K_*\) in (5i) represents the additional displacement effect produced by the wall heating.

In what follows we use the Chapman viscosity law expressed by \(\mu = C T_{31}\) for some constant C. This and additional constants such as \(\lambda , C, T_{B0}\) appearing in the above equations may be effectively removed with the aid of the following affine transformation:

$$\begin{aligned} t_*= & {} \lambda ^{-3/2}C^{-1/2}T_{B0}^{-1/2}\tau , \quad x_*=\lambda ^{-5/4} T_{B0}^{1/4}C^{-1/4} X, \\ y_*= & {} \lambda ^{-3/4}C^{1/4}T_{B0}^{-1/4} Y, \quad z_*=C^{-1/4}\lambda ^{-5/4}T_{B0}^{1/4}Z, \quad y_1=\lambda ^{-5/4}C^{-1/4}T_{B0}^{1/4} {{\bar{Y}}},\quad \\ u_{31}= & {} \lambda ^{1/4} C^{1/4}T_{B0}^{3/4} U, \quad v_*=\lambda ^{3/4} C^{3/4}T_{B0}^{1/4}V, \quad w_{31}=\lambda ^{1/4} C^{1/4}T_{B0}^{3/4} W,\\ p_{31}= & {} \lambda ^{1/2} C^{1/2}T_{B0}^{1/2} P, \quad T_{31}= T_{B0}\theta , \quad p_{11}= \lambda ^{1/2} C^{1/2}T_{B0}^{1/2} P_1,\\ A_*= & {} \lambda ^{-3/4}C^{1/4}T_{B0}^{3/4}A,\quad K_*=\lambda ^{-3/4}C^{1/4}T_{B0}^{3/4}K ,\quad D_*= C^{1/2}\lambda ^{1/2}T_{B0}^{1/2} D\\ f= & {} \lambda ^{-3/4}C^{1/4}T_{B0}^{3/4}F,\quad T_w=T_{B0}\theta _w . \end{aligned}$$

After using the above transformation in (3)–(5) in conjunction with the Chapman viscosity law and the equation of state the resulting equations are given by

$$\begin{aligned}&{\partial {U}\over \partial X}+ {\partial {V}\over \partial {Y}}+{\partial W\over \partial Z}=0 , \end{aligned}$$
(6a)
$$\begin{aligned}&{\partial {U} \over \partial \tau }+ {U}{\partial {U}\over \partial X}+ {V} {\partial {U}\over \partial {Y}}+ W{\partial U \over \partial Z}= -\theta {\partial {P}\over \partial X} + {\partial ^2 {U}\over \partial {Y}^2}, \end{aligned}$$
(6b)
$$\begin{aligned}&0= -{\partial {P}\over \partial {Y}} , \end{aligned}$$
(6c)
$$\begin{aligned}&{\partial {W} \over \partial \tau }+ {U}{\partial {W}\over \partial X}+ {V} {\partial {W}\over \partial {Y}}+ W{\partial W \over \partial Z}= -{\theta } {\partial {P}\over \partial Z} + {\partial ^2 {W}\over \partial {Y}^2}, \end{aligned}$$
(6d)
$$\begin{aligned}&{\partial {\theta } \over \partial \tau }+ {U}{\partial {\theta }\over \partial X}+ {V} {\partial {\theta }\over \partial {Y}}+W{\partial \theta \over \partial Z}= {1\over Pr}{\partial ^2 {\theta }\over \partial {Y}^2}. \end{aligned}$$
(6e)

These should be solved subject to the boundary conditions:

$$\begin{aligned}&V(\tau ,X,Y,Z)=0 \quad \mathrm{on} \quad Y=0, \end{aligned}$$
(7a)
$$\begin{aligned}&U(\tau ,X,Y,Z)=0 \quad \mathrm{on }\quad Y=0, \end{aligned}$$
(7b)
$$\begin{aligned}&W(\tau ,X,Y,Z)=0 \quad \mathrm{on }\quad Y=0, \end{aligned}$$
(7c)
$$\begin{aligned}&\theta (\tau ,X,Y,Z)= \theta _w(\tau ,X,Z) \quad \mathrm{on }\quad Y=0, \end{aligned}$$
(7d)
$$\begin{aligned}&U(\tau ,X,Y,Z)= Y \quad \mathrm{as}\quad X \rightarrow -\infty , \end{aligned}$$
(7e)
$$\begin{aligned}&U(\tau ,X,Y,Z)=Y + K(\tau ,X,Z) \quad \mathrm{as}\quad Y \rightarrow \infty , \end{aligned}$$
(7f)
$$\begin{aligned}&\theta (\tau ,X,Y,X)=1 \quad \mathrm{as} \quad Y \rightarrow \infty , \end{aligned}$$
(7g)
$$\begin{aligned}&W(\tau ,X,Y,Z) = {D \over Y} \quad \mathrm{as} \quad Y \rightarrow \infty . \end{aligned}$$
(7h)

Here K is given by

$$\begin{aligned} K(\tau ,X,Z)=F(\tau ,X,Z)+\int _0^{\infty }(\theta (\tau ,X,Y,Z)-1)\,\mathrm{d}Y + A(\tau ,X,Z) , \end{aligned}$$
(7i)

where F represents the contribution of the vibrator shape, \(\theta \) the thermal forcing and A is the displacement function. The transformed upper-deck problem is

$$\begin{aligned} (1-M_{\infty }^2){\partial ^2 P_1 \over \partial X^2}+{\partial ^2 P_1 \over \partial {{\bar{Y}}}^2}+{\partial ^2 P_1 \over \partial Z^2}=0, \end{aligned}$$

with the boundary conditions

$$\begin{aligned}&P_1 \rightarrow 0 {\quad \mathrm as \quad }(X^2+{{\bar{Y}}}^2 +Z^2) \rightarrow \infty ,\nonumber \\&P_1(\tau ,X,{{\bar{Y}}}=0,Z)=P(\tau ,X,Z),\quad {\partial P_1 \over \partial {{\bar{Y}}}}={ }{\partial ^2 A\over \partial X^2} {\quad \mathrm{on} \quad {{\bar{Y}}}=0 .} \end{aligned}$$
(8a)

In the above F represents the transformed wall shape and \(\theta _w\) is the prescribed heating profile and both these functions are assumed to be given.

The above initial-value problem is supplemented with the initial conditions

$$\begin{aligned} U=Y,\quad \theta =\theta _w=1, \quad V,W,P,P_1,A,F = 0 \quad \mathrm{for}\quad \tau \le 0 . \end{aligned}$$
(9)

The non-linear initial-value problem requires a numerical solution in general, but for small amplitudes of the vibrator and weak heating we can find a linearised solution.

2.4 Fourier–Laplace solution of the linearised equations

We will assume that the wall motion and localised heating profiles are given by

$$\begin{aligned}&F(\tau ,X,Z)= \epsilon F_a(\tau ,X,Z)=\epsilon h(X,Z)\sin (\omega _0 \tau ) , \quad \tau>0,\\&\theta _w(\tau ,X,Z)=1+\epsilon g(\tau ,X,Z), \quad \tau >0 , \end{aligned}$$

where \(\epsilon \) represents the maximum amplitude of the oscillation and \(\omega _0\) is some prescribed frequency.

For \(0<\epsilon \ll 1\) we may expand the flow quantities as

$$\begin{aligned}&U(\tau ,X,Y)= Y + \epsilon U_a(\tau ,X,Y,Z) + O(\epsilon ^2), \end{aligned}$$
(10a)
$$\begin{aligned}&V(\tau ,X,Y,Z)=\epsilon V_a(\tau ,X,Y,Z)+ O(\epsilon ^2), \end{aligned}$$
(10b)
$$\begin{aligned}&W(\tau ,X,Y,Z)=\epsilon W_a(\tau ,X,Y,Z)+ O(\epsilon ^2), \end{aligned}$$
(10c)
$$\begin{aligned}&\theta (\tau ,X,Y,Z)= 1 + \epsilon \theta _a(\tau ,X,Y,Z) + O(\epsilon ^2), \end{aligned}$$
(10d)
$$\begin{aligned}&P(\tau ,X,Z)= \epsilon P_a(\tau ,X,Z) + O(\epsilon ^2), \end{aligned}$$
(10e)
$$\begin{aligned}&P_1(\tau ,X,{{\bar{Y}}},Z)=\epsilon P_u(\tau ,X,{{\bar{Y}}},Z)+ O(\epsilon ^2), \end{aligned}$$
(10f)
$$\begin{aligned}&A(\tau ,X,Z)= \epsilon A_a(\tau ,X,Z) + O(\epsilon ^2), \end{aligned}$$
(10g)
$$\begin{aligned}&K(\tau ,X,Z)= \epsilon K_a(\tau ,X,Z) + O(\epsilon ^2). \end{aligned}$$
(10h)

Substituting (10) into (6)–(8) and linearising for small \(\epsilon \) leads to the following linearised initial-value problem:

$$\begin{aligned}&{\partial {U_a}\over \partial X}+ {\partial {V_a}\over \partial {Y}}+{\partial W_a\over \partial Z} =0 , \end{aligned}$$
(11a)
$$\begin{aligned}&{\partial {U_a} \over \partial \tau }+ {Y}{\partial {U_a}\over \partial X}+ {V_a} = -{\partial {P_a}\over \partial X} + {\partial ^2 {U_a}\over \partial {Y}^2}, \end{aligned}$$
(11b)
$$\begin{aligned}&0 = -{\partial {P_a}\over \partial {Y}} , \end{aligned}$$
(11c)
$$\begin{aligned}&{\partial {W_a} \over \partial \tau }+ {Y}{\partial {W_a}\over \partial X} =-{}{\partial P_a \over \partial Z}+ {\partial ^2 {W_a}\over \partial {Y}^2}, \end{aligned}$$
(11d)
$$\begin{aligned}&{\partial {\theta _a} \over \partial \tau }+ {Y}{\partial {\theta _a}\over \partial X} = {1\over Pr}{\partial ^2 {\theta _a}\over \partial {Y}^2} , \end{aligned}$$
(11e)
$$\begin{aligned}&K_a=h(X,Z)\sin (\omega _0 \tau )+ A_a+\int _0^{\infty }\theta _a \, \mathrm{d}Y , \end{aligned}$$
(11f)

where

$$\begin{aligned}&U_a(\tau ,X,Y=0,Z)= 0, \end{aligned}$$
(12a)
$$\begin{aligned}&V_a(\tau ,X,Y=0,Z)=0, \end{aligned}$$
(12b)
$$\begin{aligned}&W_a(\tau ,X,Y=0,Z)=0, \end{aligned}$$
(12c)
$$\begin{aligned}&\theta _a(\tau ,X,Y=0,Z)= g(\tau ,X,Z), \end{aligned}$$
(12d)
$$\begin{aligned}&U_a=V_a=W_a=\theta _a=K_a=P_a=A_a=0 \quad \mathrm{for }\quad \tau \le 0, \end{aligned}$$
(12e)
$$\begin{aligned}&U_a(\tau ,X,Y,Z)= 0 \quad \mathrm{as}\quad X \rightarrow -\infty , \end{aligned}$$
(12f)
$$\begin{aligned}&U_a(\tau ,X,Y,Z)= K_a(\tau ,X,Z) \quad \mathrm{as}\quad Y \rightarrow \infty , \end{aligned}$$
(12g)
$$\begin{aligned}&W_a(\tau ,X,Y,Z)= 0 \quad \mathrm{as}\quad Y \rightarrow \infty , \end{aligned}$$
(12h)

and

$$\begin{aligned} (1-M_{\infty }^2){\partial ^2 P_u \over \partial X^2}+{\partial ^2 P_u \over \partial {{\bar{Y}}}^2}+{\partial ^2 P_u \over \partial Z^2}=0, \end{aligned}$$
(13a)

with the boundary conditions

$$\begin{aligned}&P_u \rightarrow 0 {\quad \mathrm as \quad }(X^2+{{\bar{Y}}}^2+Z^2) \rightarrow \infty ,\nonumber \\&P_u(\tau ,X,{{\bar{Y}}}=0,Z)=P_a(\tau ,X,Z),\quad {\partial P_u \over \partial {{\bar{Y}}}}={ }{\partial ^2 A_a\over \partial X^2} {\quad \mathrm{on} \quad {{\bar{Y}}}=0 .} \end{aligned}$$
(13b)

Let us introduce the Fourier–Laplace transform

$$\begin{aligned} U_a^{\ddagger \dagger \dagger }(\omega ,k,Y,l)=\int _0^{\infty }\int _{-\infty }^{\infty } \int _{-\infty }^{\infty }U_a(\tau ,X,Y,Z) \mathrm{e}^{-\omega \tau -\mathrm{i} k X-\mathrm{i}lZ}\, \mathrm{d}X\,\mathrm{d}Z \,\mathrm{d}\tau \end{aligned}$$

and the corresponding inverse by

$$\begin{aligned} U_a(\tau ,X,Y,Z)={1\over {8\pi ^3 i}}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{{{\mathcal {L}}}} U_a^{\ddagger \dagger \dagger }(\omega ,k,Y,l)\mathrm{e}^{\omega \tau +\mathrm{i} k X+\mathrm{i}lZ}\mathrm{d}\omega \, \,\mathrm{d}l\,\mathrm{d}k , \end{aligned}$$

with similar expressions for the other quantities. The double superscript \(\dagger \dagger \) denotes the Fourier transform of both coordinates, and the superscript \(\ddagger \) the Laplace transform. The integration with respect to \(\omega \) is performed along a vertical line \({{{\mathcal {L}}}}\) in the complex \(\omega -\)plane to the right of all singularities of the transform functions to satisfy causality.

Taking transforms of (11), (12) gives

$$\begin{aligned}&\mathrm{i}k U_a^{\ddagger \dagger \dagger }+{\partial V_a^{\ddagger \dagger \dagger }\over \partial Y}+\mathrm{i}lW_a^{\ddagger \dagger \dagger }=0, \end{aligned}$$
(14a)
$$\begin{aligned}&(\mathrm{i}kY +\omega ) U_a^{\ddagger \dagger \dagger } + V_a^{\ddagger \dagger \dagger }=-\mathrm{i}kP_a^{\ddagger \dagger \dagger }+ {\partial ^2 U_a^{\ddagger \dagger \dagger } \over \partial Y^2} . \end{aligned}$$
(14b)
$$\begin{aligned}&(\mathrm{i}kY +\omega ) W_a^{\ddagger \dagger \dagger } =-{\mathrm{i}l}P_a^{\ddagger \dagger \dagger }+ {\partial ^2 W_a^{\ddagger \dagger \dagger } \over \partial Y^2} . \end{aligned}$$
(14c)
$$\begin{aligned}&(\mathrm{i}kY +\omega ) \theta _a^{\ddagger \dagger \dagger } ={1\over Pr} {\partial ^2 \theta _a^{\ddagger \dagger \dagger } \over \partial Y^2} . \end{aligned}$$
(14d)

Equation (14d) for the temperature perturbation may be solved in terms of Airy functions to obtain the solution

$$\begin{aligned} \theta _a^{\ddagger \dagger \dagger }= D_0\, \mathrm{Ai}(Pr^{1/3}\xi )+D_1\,\mathrm{Bi}(Pr^{1/3}\xi ) , \end{aligned}$$
(15)

where

$$\begin{aligned} \xi =(\mathrm{i}k)^{1/3}Y + \xi _0, \quad \xi _0=\omega (\mathrm{i}k)^{-{2/3}} . \end{aligned}$$
(16)

We will take a branch cut along the positive imaginary axis in the complex plane so that \(-3\pi /2<\arg (k)< \pi /2\). Then the function \(\mathrm Bi(\xi )\) grows exponentially when \(Y \rightarrow \infty \) and hence \(D_1\) must be zero. Application of the boundary conditions yields

$$\begin{aligned} \theta _a^{\ddagger \dagger \dagger }= g^{\ddagger \dagger \dagger }(\omega ,k,l){\mathrm{Ai}(Pr^{1/3}\xi ) \over \mathrm{Ai}(\eta _0)}, \end{aligned}$$
(17)

where we have written

$$\begin{aligned} \eta _0=Pr^{1/3}\xi _0. \end{aligned}$$
(18)

Next introduce

$$\begin{aligned} Q_a^{\ddagger \dagger \dagger }=k U_a^{\ddagger \dagger \dagger }+l W_a^{\ddagger \dagger \dagger }, \end{aligned}$$

so that from (14b), (14c) we obtain

$$\begin{aligned} (\mathrm{i}kY+\omega )Q_a^{\ddagger \dagger \dagger }+ kV_a= -\mathrm{i}(k^2+l^2)P_a+ {\partial ^2 Q_a^{\ddagger \dagger \dagger } \over \partial Y^2}. \end{aligned}$$
(19)

Differentiating (19) with respect to Y and using the continuity equation gives

$$\begin{aligned} {\partial ^3 Q_a^{\ddagger \dagger \dagger }\over \partial Y^3}-(\mathrm{i}kY+\omega ){\partial Q_a^{\ddagger \dagger \dagger } \over \partial Y}=0 . \end{aligned}$$

This has the solution

$$\begin{aligned} {\partial Q_a^{\ddagger \dagger \dagger }\over \partial Y}= C_0\, \mathrm{Ai}(\xi )+ C_1\, \mathrm{Bi}(\xi ). \end{aligned}$$
(20)

The Airy function \(\mathrm Bi(\xi )\) grows exponentially for large Y and so by the same reasoning as earlier we must take \(C_1=0\). Setting \(Y=0\) in (19) and using (20) gives

$$\begin{aligned} (\mathrm{i}k)^{1/3}C_0 \mathrm{Ai}'(\xi _0)= \mathrm{i}(k^2+{l^2}) P_a^{\ddagger \dagger \dagger } . \end{aligned}$$
(21)

We can further integrate (20) to obtain

$$\begin{aligned} Q_a^{\ddagger \dagger \dagger }=C_0(\mathrm{i}k)^{-{1\over 3}}\int _{\xi _0}^{\xi }\mathrm{Ai}(\xi )\,\mathrm{d}\xi . \end{aligned}$$
(22)

Letting \(Y\rightarrow \infty \) in (22) and using the transformed boundary conditions from (12) shows that

$$\begin{aligned} kK_a^{\ddagger \dagger \dagger }= C_0 (ik)^{-{1\over 3}}\int _{\xi _0}^{\infty }\mathrm{Ai}(\xi )\,\mathrm{d}\xi . \end{aligned}$$
(23)

We also have from (11f)

$$\begin{aligned} K_a^{\ddagger \dagger \dagger }= A_a^{\ddagger \dagger \dagger }+\int _{0}^{\infty }\theta _a^{\ddagger \dagger \dagger }\,\mathrm{d}Y +{h^{\dagger \dagger }(k,l) \omega _0 \over (\omega ^2+\omega _0^2)}, \end{aligned}$$
(24)

which after using the solution for \(\theta _a^{\ddagger \dagger \dagger }\) becomes

$$\begin{aligned} K_a^{\ddagger \dagger \dagger }= A_a^{\ddagger \dagger \dagger }+g^{\ddagger \dagger \dagger }(\mathrm{i}kPr)^{-1/3}\int _{\eta _0}^{\infty }{\mathrm{Ai}(\eta ) \over \mathrm{Ai(\eta _0)}}\,\mathrm{d}\eta +{h^{\dagger \dagger }(k,l) \omega _0 \over (\omega ^2+\omega _0^2)} . \end{aligned}$$
(25)

Equations (13a)–(13b) for \(P_u\) do not involve \(\tau \) explicitly and therefore taking Fourier–Laplace transforms of (13a) and applying the boundary conditions gives the usual relation,

$$\begin{aligned} P_a^{\ddagger \dagger \dagger }={k^2 \over (k^2(1-M_{\infty }^2)+l^2)^{1\over 2}}A_a^{\ddagger \dagger \dagger } . \end{aligned}$$
(26)

Finally eliminating \(C_0\) and solving for \(P_a^{\ddagger \dagger \dagger }\) from (21), (23), (26) gives

$$\begin{aligned} P_a^{\ddagger \dagger \dagger }(\omega ,k,l)=P_{\pm }^{\ddagger \dagger \dagger }(\omega ,k,l) = {H^{\ddagger \dagger \dagger }(\omega ,k,l) \omega _0 |k| \mathrm{Ai}'(\xi _0) \over (\omega ^2+\omega _0^2) D^{\pm }(\xi _0,k)}, \end{aligned}$$
(27)

where \(\xi _0\) is defined in (16) and

$$\begin{aligned} D^{\pm }(\xi _0,k)= -{(k^2\beta ^2+l^2)^{1\over 2}\over |k|}\mathrm{Ai}'(\xi _0) \mp k(ik)^{-{5\over 3}} (k^2+{l^2 })\int _{\xi _0}^{\infty }\mathrm{Ai}(\xi )\,\mathrm{d}\xi . \end{aligned}$$
(28)

Here we have put \(\beta ^2 = 1-M_{\infty }^2\) and the plus sign in \(D^{\pm }\) corresponds to k positive and the minus sign to k negative. In the above expression we have defined

$$\begin{aligned} H^{\ddagger \dagger \dagger }(\omega ,k,l)=h^{\dagger \dagger }(k,l)+ g^{\ddagger \dagger \dagger }(\omega ,k,l){(\omega ^2+\omega _0^2)(\mathrm{i}k Pr)^{-1/3}\int _{\eta _0}^{\infty } \mathrm{Ai}(\eta )\,\mathrm{d}\eta \over \omega _0 \mathrm{Ai}(\eta _0)}. \end{aligned}$$
(29)

The disturbed pressure \(P_a(\tau ,X)\) is calculated by formally inverting (27).

A number of results are immediately apparent from (27). Notice that if we have no localised heating and set \(g^{\ddagger \dagger \dagger }\) to be zero, then (27) reduces to the expression obtained by [14, 16] for the two-dimensional case and [20] for the three-dimensional problem.

We notice from (27) that even without a vibrator, localised heating can excite Tollmien–Schlichting disturbances. With a vibrator present, if the localised heating profile is chosen such that \(H^{\ddagger \dagger \dagger }=0\), then the response \(P_a^{\ddagger \dagger \dagger }\) is zero, which means no Tollmien–Schlichting disturbances. In fact the required localised heating profile is given by \(g^{\ddagger \dagger \dagger }(\omega ,k)=g_{TC}^{\ddagger \dagger \dagger }\) where

$$\begin{aligned} g_{TC}^{\ddagger \dagger \dagger }(\omega ,k,l)=-{h^{\dagger \dagger }(k,l) \omega _0 \,\mathrm{Ai}(\eta _0)(\mathrm{i}k Pr)^{1/3} \over (\omega ^2+\omega _0^2)\int _{\eta _0}^{\infty }\mathrm{Ai}(\eta )\,\mathrm{d}\eta } . \end{aligned}$$
(30)

The expression (30) is an important new result showing that it is possible to control Tollimien–Schlichting instabilities in the boundary layer. It is the three-dimensional counterpart to the result obtained by [14] with \(h^{\dagger \dagger }(k,l)\) replacing the \(h^{\dagger }(k)\) in their problem. In the numerical solutions presented below, an approximation to the cancellation function (30) is used successfully to substantially reduce the amplitude of unstable Tollmien–Schlichting disturbances generated by the vibrator.

3 Numerical solution of the linearised triple-deck initial-value problem

Whilst it is possible to obtain an analytical solution of the initial-value problem, it is more convenient to obtain results using a numerical solution.

In the numerical work we choose the wall-shape and heating functions to be

$$\begin{aligned} F_a(\tau ,X,Z)=h(X,Z) q(\tau ) , \quad \theta _{aw}=g(\tau ,X,Z) \quad \mathrm{for} \quad \tau >0 . \end{aligned}$$
(31)

The function h(XZ) is a smooth Gaussian hump given by

$$\begin{aligned} h(X,Z)= \mathrm{e}^{-{\pi } (\delta _x^2X^2+ \delta _z^2Z^2)}, \end{aligned}$$

and for \(q(\tau )\) we have taken

$$\begin{aligned} q(\tau )=(1-\mathrm{e}^{-a\tau ^2})\sin (\omega _0 \tau ). \end{aligned}$$
(32)

In the numerical computations presented below we have taken \(\delta _x=\delta _z=1/4\) and for the function \(q(\tau )\), \(a=1/10\) which gives a smoother initial start. The choice of the wall temperature profile g is discussed below.

Equation (11) were solved in Fourier transform space with a second-order time marching scheme. The details of the numerical techniques used are very similar to those used in the two-dimensional problem looked at by [14] and therefore only brief details and relevant changes are highlighted.

First we can write Eqs. (11)–(13) in transform space as

$$\begin{aligned}&\mathrm{i}kU_a^{\dagger \dagger } + {\partial V_a^{\dagger \dagger } \over \partial Y} +\mathrm{i} l W_a^{\dagger \dagger }=0, \end{aligned}$$
(33a)
$$\begin{aligned}&{\partial U_a^{\dagger \dagger } \over \partial \tau } + \mathrm{i}kY U_a^{\dagger \dagger } + V_a^{\dagger \dagger } = -\mathrm{i}k P_a^{\dagger \dagger } + {\partial ^2 U_a^{\dagger \dagger } \over \partial Y^2}, \end{aligned}$$
(33b)
$$\begin{aligned}&{\partial W_a^{\dagger \dagger } \over \partial \tau } + \mathrm{i}kY W_a^{\dagger \dagger } = -\mathrm{i}l P_a^{\dagger \dagger } + {\partial ^2 W_a^{\dagger \dagger } \over \partial Y^2}, \end{aligned}$$
(33c)
$$\begin{aligned}&{\partial \theta _a^{\dagger \dagger } \over \partial \tau } + \mathrm{i}kY \theta _a^{\dagger \dagger } = {1\over Pr}{\partial ^2 \theta _a^{\dagger \dagger } \over \partial Y^2}, \end{aligned}$$
(33d)
$$\begin{aligned}&U_a^{\dagger \dagger }(\tau ,k,Y=0,l)=V_a(\tau ,k,Y=0,l)=W_a(\tau ,k,Y=0,l)=0, \end{aligned}$$
(33e)
$$\begin{aligned}&\theta _a(\tau ,k,Y=0,l)=g^{\dagger \dagger }(\tau ,k,l), \end{aligned}$$
(33f)
$$\begin{aligned}&K_a^{\dagger \dagger }(\tau ,k,l)=q(\tau )h^{\dagger \dagger }(k,l)+A_a^{\dagger \dagger } (\tau ,k,l)+\int _0^{\infty }\theta _a^{\dagger \dagger }(\tau ,k,Y,l)\,\mathrm{d}Y , \, \tau >0 . \end{aligned}$$
(33g)
$$\begin{aligned}&U_a^{\dagger \dagger }(\tau ,k, Y\rightarrow \infty ,l)= K_a^{\dagger \dagger }(\tau ,k,l), \end{aligned}$$
(33h)
$$\begin{aligned}&P_a^{\dagger \dagger }={k^2 \over (k^2(1-M_{\infty }^2)+l^2)^{1\over 2}}A_a^{\dagger \dagger } . \end{aligned}$$
(33i)

By introducing the cross-flow variable

$$\begin{aligned} Q_a^{\dagger \dagger }(\tau ,k,Y,l)= k U_a^{\dagger \dagger }(\tau ,k,Y,l)+ l V_a^{\dagger \dagger }(\tau ,k,Y,l) \end{aligned}$$

and using a second-order fully implicit time-differencing scheme with time step \(d\tau \), for each k and l the system of equations in Y are solved using a Chebychev collocation method. This leads to a linear system of the form:

$$\begin{aligned} \mathbf{G}\left( \begin{array}{c}\mathbf{Q}_a^{{\dagger \dagger }(n+1)}\\ A_a^{{\dagger \dagger }(n+1)}\\ {\varvec{\theta }}_a^{{\dagger \dagger }(n+1)}\end{array}\right) =\mathbf{R}^{(n)} \end{aligned}$$
(34)

for the unknowns \(\mathbf{Q}_a^{\dagger \dagger (n+1)},A_a^{\dagger \dagger (n+1)},{\varvec{\theta }}_a^{\dagger \dagger (n+1)}.\) Here the elements

$$\begin{aligned} Q_{aj}^{\dagger \dagger (n)}=Q_a^{\dagger \dagger }(\tau _n, k,Y_j,l) \end{aligned}$$

of the vector \(\mathbf{Q}_a^{\dagger \dagger (n)}\) are defined at the collocation points

$$\begin{aligned} Y=Y_j= {Y_{max} \over 2}\left( 1-\cos \left( {j \pi \over N}\right) \right) , \quad j=0, \ldots , N , \end{aligned}$$

with similar expressions for the other quantities, and \(\tau _n=ndt\). The approximate outer boundary \(Y_{max}=40\) was typically used in our computations. For a given kl pair we can solve the linear system (34) to find the unknown variables and in particular \( P_a^{\dagger \dagger (n)}(k,l)\). This is then inverted using the discrete inverse Fast Fourier Transform to obtain the pressure \(P_a(\tau _n,X,Z)\) as well as the velocity, temperature and displacement functions. For the results shown below we have used 1024 k modes, 128 modes in l, with 256 points per time period \(T_{p}=2\pi /\omega _0\), and with 32 Chebychev modes. We have also used \(Pr=1\). Extensive grid size and other checks have been carried out. Previous work [14] mentions the difficulties which arise from using too large a timestep and from the choice of the range of wave numbers used in the computations.

4 Results

The dispersion relation for the three-dimensional problem is independent of forcing or thermal effects and has been analysed by [19, 20]. The salient properties are that for \(\omega \) real and positive, all the roots are situated in the left-hand side of the complex \(k-\)plane and there are countably many roots of which only one is unstable.

The dispersion relation for the 3D problem is given by setting \(D_+=0\) for k positive and \(D_{-} = 0\) for k negative in Eq. (28). As in [19, 20], we reduce this relation to its 2D analogue by introducing the variable

$$\begin{aligned} k' = \text {sign}(k) |k|^{\frac{1}{4}} (k^2+l^2)^{\frac{3}{4}}(k^2\beta ^2+l^2)^{-\frac{3}{8}} . \end{aligned}$$
(35)

As a result (35) assumes the form which has been studied for example by [16, 26], that is

$$\begin{aligned} \frac{\text {Ai}'(\xi _0')}{\int _{\xi _0'}^{\infty } \text {Ai}(\xi )\mathrm{d} \xi } \pm |k'|^2 k' (\mathrm{i}k')^{-\frac{5}{3}}=0, \end{aligned}$$
(36)

with \(\xi _0' = \omega ' (\mathrm{i}k')^{-\frac{2}{3}}\). The dispersion relation has roots given by \(\omega _n' = (\mathrm{i}k')^{\frac{2}{3}} \xi _{0n}'\) with \(\omega _1'\) being the only unstable root. For \(|k'| > k'^*\), \(\mathfrak {R}(\omega _1')>0\) indicating instability in time. Similarly for complex wavenumbers, past this point \(\mathfrak {I}(k')<0\) which also shows instability. Here the value of the critical point \(k'=k'^*\) is dependant upon the Prandtl–Glauert factor \(\beta \), and thus the Mach number \(M_\infty \) [26]. Savenkov [20] introduced a generalisation of the Squire transformation,

$$\begin{aligned} \omega _1(k,l) = \Bigg ( \frac{k}{k'}\Bigg )^{\frac{2}{3}}\omega _1'(k'), \end{aligned}$$
(37)

which allows us to study the influence of changing \(M_\infty \), and thus the compressibility of the flow, on the flow stability. Moreover, by introducing the variable \({\hat{\beta }} = l/k\), [20] showed that we can find \(\max (\mathfrak {R}(\omega _1(k,l)))\), the maximum growth rate, by writing

$$\begin{aligned}&\omega _1(k,l) = f({\hat{\beta }})\omega _1'(k'), \end{aligned}$$
(38)
$$\begin{aligned}&f({\hat{\beta }}) = (1+{\hat{\beta }}^2)^{-\frac{1}{2}}({\hat{\beta }}^2 + (1-M_\infty ^2))^{\frac{1}{4}}. \end{aligned}$$
(39)

Now

$$\begin{aligned} \max (\mathfrak {R}(\omega _1)) = f_e \sigma _{0e}, \end{aligned}$$
(40)

where \(f_e = \max (f({\hat{\beta }}))\) and \(\sigma _{0e} = \max \mathfrak {R}(\omega _1')\) which is achieved at \(|k'| = k_2'^*\) and, like \(k'^*\) is dependant upon the Prandtl–Glauert factor. It can be shown that \(f_e\) occurs at \({\hat{\beta }} = {\hat{\beta }}_e\) with \(k = k_e = k_2'^*f_e^{\frac{3}{2}}\) and \(l=l_e = {\hat{\beta }}_e k_e\). In particular the important results due to [20] are

$$\begin{aligned}&{\hat{\beta }}_e = 0,\quad k_e = (1-M_\infty ^2)^{\frac{3}{8}} k_2'^*,\quad l_e = 0 \quad \text {for}\quad M_\infty \le M_\infty ^*, \end{aligned}$$
(41)
$$\begin{aligned}&{\hat{\beta }}_e = (2M_\infty ^2 - 1)^{\frac{1}{2}},\quad k_e = (2 M_\infty )^{-\frac{3}{4}} k_2'^*, \nonumber \\&l_e = (2 M_\infty ^2 - 1)^{\frac{1}{2}} k_e \quad \text {for}\quad M_\infty \ge M_\infty ^*, \end{aligned}$$
(42)

where \(M_\infty ^* = 1/\sqrt{2}\).

Figure 2 shows how \(\max (\mathfrak {R}(\omega _1))\) varies with \(M_\infty \). In Fig. 3 we see that keeping l constant and increasing \(M_\infty \) lowers \(\max (\mathfrak {R}(\omega _1))\). Also, we see in Fig. 4 that as \(l\rightarrow l_e\) from above or below, the maximum growth rate \(\max (\mathfrak {R}(\omega _1))\) increases. It can be seen that for k positive, the results shown in Fig. 3b can be obtained via a mirror reflection of the results for k negative in Fig. 3a. The same observation also applies to Fig. 4. In Fig. 5 we have plotted contours of \(\mathfrak {R}(\omega _1(k,l)/\max (\mathfrak {R}(\omega _1(k,l))\) for \(M_{\infty }=0.5\) and \(M_{\infty }=0.85\). The results for \(M_{\infty }=0.85\) are in excellent agreement with the corresponding figure in [20]. The neutral curve is shown by the zero contour. The maximum growth rate for \(M_{\infty }=0.85\) is for an oblique mode and shown at the point labelled 1.

In summary, if \(M_\infty < M_\infty ^*=1/\sqrt{2}\) the two-dimensional mode with \(l=0\) is the most unstable and instability arises at the critical frequency \(\omega _0\) with corresponding neutral wave-number \(k = k_1^*\) which is dependant upon the Mach number. For \(M_\infty > M_\infty ^*\), the oblique mode (\(l\ne 0\)) is most unstable.

Fig. 2
figure 2

For \(l=l_e\), \(k=k_e\), the maximum growth rate \(\max (\mathfrak {R}(\omega _1)) = f_e \sigma _{0e}\) varies with \(M_\infty \) according to how \(f_e\) and \(\sigma _{0e}\) vary with \(M_\infty \). Note that for \(M_{\infty }<M_{\infty }^*\), \(l_e=0\)

Fig. 3
figure 3

The unstable root \(\omega _1\) of \(D_{-}=0\) (a) and \(D_{+}=0\) (b) for different values of \(M_\infty \) with \(l =0.4\) fixed and k varying from \((- \infty ,0)\) (a) and \((0,\infty )\) (b)

Fig. 4
figure 4

The unstable root \(\omega _1\) of \(D_{-}=0\) (a) and \(D_{+}=0\) (b) for different values of l with \(M_\infty =0.75\) fixed and k varying from \((- \infty ,0)\) (a) and \((0,\infty )\) (b)

Fig. 5
figure 5

Contours of \(\mathfrak {R}(\omega _1)/\max (\mathfrak {R}(\omega _1))\) as a function of kl for a \(M_{\infty }=0.5\) and b \(M_{\infty }=0.85\)

The numerical solution of the initial-value problem for \(M_{\infty }=0\) confirms the analytical predictions. In Figs. 6, 7, 8 we show the development of the pressure response at integer multiples of the time period for the non-heated case at a stable frequency \(\omega =2\) in Fig. 6, neutral frequency \(\omega =2.298\) in Fig. 7 and unstable frequency \(\omega =2.5\) in Fig. 8.

In Figs. 9 and  10 we show similar results taking \(M_{\infty }=0.75\) and \(M_{\infty }=0.85\) with a forcing frequency \(\omega =2.5\). As compared to the incompressible case, non-zero Mach numbers reduces the lateral development of the disturbance but with the Mach number approaching \(M_{\infty }=1^-\) the disturbance growth in the oblique direction is clearly visible in Fig. 10. Also the disturbance amplitudes for the non-zero Mach numbers are also much larger as compared to the incompressible case.

Fig. 6
figure 6

Pressure response \(P_a(\tau ,X,Z)\) at times \(\tau =2T_p,4T_p,6T_p, 8T_p\) for the stable case with no heating. Here \(T_p\) is the period of oscillation of the vibrator, \(\omega =2\) and \(M_{\infty }=0\)

Fig. 7
figure 7

Pressure response \(P_a(\tau ,X,Z)\) at times \(\tau =2T_p,4T_p,6T_p, 8T_p\) for the neutral case with \(\omega =2.298\) and no heating. Here \(T_p\) is the period of oscillation of the vibrator and \(M_{\infty }=0.\)

Fig. 8
figure 8

Pressure response \(P_a(\tau ,X,Z)\) at times \(\tau =2T_p,4T_p,6T_p, 8T_p\) for the unstable case with \(\omega =2.5\) and no heating. \(T_p\) is the period of oscillation of the vibrator and \(M_{\infty }=0\)

Fig. 9
figure 9

Pressure response \(P_a(\tau ,X,Z)\) at times \(\tau =2T_p,4T_p,6T_p, 8T_p\) for the unstable case with \(\omega =2.5\) with \(M_{\infty }=0.75\) and no heating. Here \(T_p\) is the period of oscillation of the vibrator

Fig. 10
figure 10

Pressure response \(P_a(\tau ,X,Z)\) at times \(\tau =2T_p,4T_p,6T_p, 8T_p\) for the unstable case with \(\omega =2.5\) with \(M_{\infty }=0.85\) and no heating. Here \(T_p\) is the period of oscillation of the vibrator

Fig. 11
figure 11

A comparison of the pressure response for the unstable case with \(\omega =2.5\) and \(M_{\infty }=0\) at integer multiples of the time period \(T_p\). On the left is the response without heating, and the right-hand figure shows the same with the cancellation formula (43) applied. a At \(\tau =3T_p\), b \(\tau =5T_p\), c \(\tau =7T_p\)

Fig. 12
figure 12

A comparison of the pressure response for the unstable case with \(\omega =2.5\) and \(M_{\infty }=0.85\) at integer multiples of the time period \(T_p\). On the left is the response without heating, and the right-hand figure shows the same with the cancellation formula applied. a At \(\tau =3T_p\), b \(\tau =5T_p\), c \(\tau =7T_p\)

Fig. 13
figure 13

Contour plots of the vibrator motion (left side figures) and the surface heating control according to the formula (43) applied (right side figures) at various time intervals. a At \(\tau =7T_p\), b \(\tau =7{1\over 4}T_p\), c \(\tau =7{1\over 2}T_p\), d \(\tau =7{3\over 4} T_p\)

4.1 Tollmien–Schlichting wave cancellation

The result (30) for the heating function \(g_{TC}\) required to cancel the Tollmien–Schlichting waves can be used to derive an approximate formula which we have confirmed in our computations to significantly reduce the amplitude of the Tollmien–Schlichting disturbances generated by the vibrator. The inversion for the two-dimensional case is discussed in [14] and apart from the change in the shape of the hump represented by \(h^{\dagger \dagger }\) there are few differences. In fact the approximation \(g_{TCN}(\tau ,X,Z)\) used in our numerical work is given by

$$\begin{aligned} g^{{\dagger \dagger }}_{TCN}(\tau ,k,l)= & {} {h^{{\dagger \dagger }}(k,l) {(\mathrm{i}kPr)^{1/3}\mathrm{Ai}(\eta _0(-\mathrm{i}\omega _0,k)) \over 2\mathrm{i} \int _{\eta _0(-\mathrm{i}\omega _0,k)}^{\infty } \mathrm{Ai}(\eta )\,\mathrm{d}\eta } }\mathrm{e}^{-\mathrm{i}\omega _0\tau } -{h^{{\dagger \dagger }}(k,l) {\mathrm{Ai}(\eta _0(\mathrm{i}\omega _0,k))(\mathrm{i}kPr)^{1/3} \over 2\mathrm{i} \int _{\eta _0(\mathrm{i}\omega _0,k)}^{\infty } \mathrm{Ai}(\eta )\,\mathrm{d}\eta }}\mathrm{e}^{\mathrm{i}\omega _0\tau } . \end{aligned}$$
(43)

We note that the approximate formula (43) includes the dominant contributions from the poles of denominator in (30) at the forcing frequency \(\omega _0\). The additional contributions stemming from the (countable) zeros of the generalised Airy function \(\int _{\eta _0}^{\infty } \mathrm{Ai}(\eta ) \, \mathrm{d}\eta \) in (30) have terms which are exponentially decreasing in time.

In Fig. 11 we show this formula in action for an unstable frequency of \(\omega =2.5\) at various instances in time. The heating profile used was the function \(g_{TCN}\). Note that even though in the analysis we used the wall motion function \(F_a(\tau ,X,Z)=h(X,Z)\sin (\omega _0\tau )\), the numerical computations are performed with the smooth initial start given by (32). After a short time, the wave cancellation takes effect and as Fig. 11c shows the Tollmien–Schlichting wave amplitude is significantly diminished. Surface heating on its own is able to generate Tollmien–Schlichting waves, and indeed can reinforce the disturbance as Fig. 11a shows, but with an appropriately chosen wall heating profile the wave generated by the vibrating hump shape can be eliminated.

Even though the formula (43) does not have explicit dependence on the Mach number, it also works for the compressible case. In Fig. 12 we show a similar comparison for \(M_{\infty }=0.85\) between the results for Tollmien–Schlichting wave generation by the vibrator with and without the cancellation formula. Again, as seen clearly in Fig. 12b, c after a short time the cancellation is effective in substantially reducing the unstable wave amplitude. In fact, although not shown here, the formula is equally effective also for the stable and neutral cases.

In Fig. 13 we have compared the wall motion with the contour plots of the heating function used based on the approximate formula (43) at various time intervals. The heating/cooling is localised and mostly out of phase with the vibrator motion, with a certain degree of overcompensation. As Fig. 13a, c show, at instants when the vibrator motion is zero, the heating function is clearly not zero.

5 Conclusions

Control of unstable disturbances in boundary layer flows is of considerable interest especially with regard to delaying transition. In quiet disturbance environments the Tollmien–Schlichting instability has a major role in inducing early transition. In this paper we have shown how it is possible to introduce unsteady surface heating to substantially reduce Tollmien–Schlichting disturbance amplitudes. The formula (43) is a generalisation of the corresponding formula for two-dimensional flows obtained by Brennan et al. [14] and our numerical simulations have shown how effective this is for a variety of different compressible flows. The formula was derived based on linear theory and it would be useful to develop this further by implementing the same ideas but for the full Navier–Stokes equations. The extension of the current work to fully non-linear simulations of the unsteady triple-deck equations would be of considerable interest but is beyond the scope of this paper. In principle the numerical technique used here could be modified to handle the non-linear terms as in the work of [27] for instance, by evaluating the omitted non-linear terms in physical space and transforming back to Fourier space. For the triple-deck type equations, this has been done for predominantly steady-state flows but there has been limited application of the ideas to perform unsteady simulations.

Again, as in [14], the cancellation formula we have derived is somewhat special in that it applies to a vibrator placed at the wall and as noted by a referee, the generation and suppression of the disturbances is fully managed. In reality disturbances can be generated via many different sources and it would be interesting to be able to obtain a corresponding result for more general disturbance sources. This is beyond the scope of the current paper, but in principle the basic ideas could be incorporated in the analysis of receptivity. This requires further investigation.

We have used an asymptotic approach to modelling three-dimensional disturbance growth in boundary layers. Whilst it is not possible to make direct comparisons with the experimental work of [7, 12], some qualitative comparisons between our results and those of Gaster [28] can be made. For example Figs. 9, 10 of [28] shows experimental results for a wave-packet travelling downstream and generated by a pulsed disturbance source on a flat plate. This can be compared with Figs. 8 and 11 in our work. The wave-packet shapes look very similar and provides encouragement for utilising our results for the compressible boundary layers which are much harder to study experimentally.