2.1 Scalings and governing equations

We consider the boundary-layer flow over a concave wall with a characteristic radius of curvature $r_{0}^{\ast }$ . The incompressible oncoming flow is assumed to be uniform with a speed $U_{\infty }$ , superimposed on which are small-amplitude vortical fluctuations. The FSVD has a characteristic spanwise length scale $\unicode[STIX]{x1D6EC}$ . The analysis is pertinent to low-frequency, including zero-frequency components, which can entrain into the boundary layer to excite Görtler vortices. For simplicity, we consider the case of a pair of free-stream vortical modes with opposite (dimensional) spanwise wavenumbers $\pm k_{3}^{\ast }$ . Although this is a drastic idealization of physical reality, the assumed form of disturbances suffices for elucidating key physical mechanisms and transition process. The generalization of the formulation and calculation to realistic broadband FSVD is rather straightforward.

We describe the flow in a curvilinear system $(x^{\ast },y^{\ast },z^{\ast })$ with its origin at the leading edge, where $x^{\ast }$ measures the distance to the leading edge along the wall, $y^{\ast }$ is the distance to the wall and $z^{\ast }$ the spanwise coordinate normal to both $x^{\ast }$ and $y^{\ast }$ . Let

(2.1a,b ) $$\begin{eqnarray}\left(x,y,z\right)=(x^{\ast },y^{\ast },z^{\ast })/\unicode[STIX]{x1D6EC},\quad t=U_{\infty }t^{\ast }/\unicode[STIX]{x1D6EC},\end{eqnarray}$$

where $t^{\ast }$ is the dimensional time. The Reynolds number is defined as

(2.2) $$\begin{eqnarray}R_{\unicode[STIX]{x1D6EC}}=U_{\infty }\unicode[STIX]{x1D6EC}/\unicode[STIX]{x1D708}\end{eqnarray}$$

and the global Görtler number as (WZL)

(2.3) $$\begin{eqnarray}G_{\unicode[STIX]{x1D6EC}}=R_{\unicode[STIX]{x1D6EC}}^{2}\unicode[STIX]{x1D6EC}/r_{0}^{\ast }.\end{eqnarray}$$

We assume that $R_{\unicode[STIX]{x1D6EC}}$ is asymptotically large, i.e. $R_{\unicode[STIX]{x1D6EC}}\gg 1$ , but take $G_{\unicode[STIX]{x1D6EC}}$ as being of $O(1)$ so that the resulting formulation accommodates all regimes of Görtler instability.

Among the Fourier components of FSVD, of importance are those with low frequencies or equivalently long streamwise wavelengths $2\unicode[STIX]{x03C0}/k_{1}^{\ast }\gg \unicode[STIX]{x1D6EC}$ , where $k_{1}^{\ast }$ is the dimensional wavenumber. Specifically, when $2\unicode[STIX]{x03C0}/k_{1}^{\ast }=O(\unicode[STIX]{x1D6EC}R_{\unicode[STIX]{x1D6EC}})$ , the induced response in the boundary layer takes on the character of Görtler vortices at distances $x^{\ast }\sim \unicode[STIX]{x1D6EC}R_{\unicode[STIX]{x1D6EC}}$ , as was shown by WZL. We thus introduce the slow streamwise and time variables,

(2.4a,b ) $$\begin{eqnarray}\hat{x}=x/R_{\unicode[STIX]{x1D6EC}},\quad \hat{\unicode[STIX]{x1D70F}}=t/R_{\unicode[STIX]{x1D6EC}}.\end{eqnarray}$$

The disturbances in the upstream region are passively advected by the uniform background flow $U_{\infty }$ . It follows that the total velocity field can be written as $(u^{\ast },v^{\ast },w^{\ast })/U_{\infty }=(1,0,0)+\tilde{\boldsymbol{u}}_{\infty }$ , and the disturbance $\tilde{\boldsymbol{u}}_{\infty }$ takes the form

(2.5) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{\infty }=\unicode[STIX]{x1D716}\boldsymbol{u}_{\infty }(\hat{x}-\hat{\unicode[STIX]{x1D70F}},y,z)=\unicode[STIX]{x1D716}(\hat{\boldsymbol{u}}_{+}^{\infty }\text{e}^{\text{i}k_{3}z}+\hat{\boldsymbol{u}}_{-}^{\infty }\text{e}^{-\text{i}k_{3}z})\text{exp}[\text{i}\hat{k}_{1}(\hat{x}-\hat{\unicode[STIX]{x1D70F}})+\text{i}k_{2}y]+\text{c.c}.,\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is a measure of the disturbance intensity, $\hat{\boldsymbol{u}}_{\pm }^{\infty }=\{\hat{u} _{1,\pm }^{\infty },\hat{u} _{2,\pm }^{\infty },\hat{u} _{3,\pm }^{\infty }\}=O(1)$ is the scaled velocity, $\boldsymbol{k}=\{k_{1},k_{2},k_{3}\}$ is the wavenumber vector and $\hat{k}_{1}=k_{1}R_{\unicode[STIX]{x1D6EC}}$ . The continuity condition implies that

(2.6) $$\begin{eqnarray}k_{1}\hat{u} _{1,\pm }^{\infty }+k_{2}\hat{u} _{2,\pm }^{\infty }\pm k_{3}\hat{u} _{3,\pm }^{\infty }=0.\end{eqnarray}$$

The turbulent Reynolds number is defined as (Goldstein et al. Reference Goldstein, Leib and Cowley1992, WZL)

(2.7) $$\begin{eqnarray}r_{t}\equiv \unicode[STIX]{x1D716}R_{\unicode[STIX]{x1D6EC}}=O(1).\end{eqnarray}$$

It is necessary here to make a few remarks on FSVD and the so-called ‘free-stream coherent structures’, which have been described in recent studies; see for example Deguchi & Hall (Reference Deguchi and Hall2017) and references therein. The former are external disturbances characterizing the flow environment and present in the main inviscid part of the flow. In contrast, the latter are the nonlinearly self-sustained intrinsic disturbances, and decay to zero in the main inviscid region although they reside in the outer reach of the boundary layer, in the so-called ‘production layer’. FSVD and free-stream coherent structures are therefore different physical entities. Mathematically, there is a certain degree of similarity: both have a phase speed close to the free-stream velocity, and short-wavelength FSVD have an ‘edge layer’ (Dong & Wu Reference Dong and Wu2013; Wu & Dong Reference Wu and Dong2016), which is essentially the same as the ‘production layer’ of free-stream coherent structures. Recently, Deguchi & Hall (Reference Deguchi and Hall2017) showed that the nonlinear interactions of free-stream coherent structures in the production layer generate streaks, which can in turn excite Görtler vortices in the boundary layer if the wall is concave. We note that a similar excitation can take place for short-wavelength FSVD too. Presently, we consider long-wavelength FSVD, which are a more relevant and efficient agent for generating Görtler vortices.

Figure 1. Schematic illustration of the physical problem and the asymptotic structure.

Similar to the flat-plate case considered by LWG, the ensuing boundary layer over a concave wall evolves through four asymptotic regimes or regions (WZL), illustrated in figure 1. Near the leading edge is an inviscid region I, which is $O(\unicode[STIX]{x1D6EC})$ in all three directions, and the disturbance is governed by the linear rapid-distortion theory (LWG). Beneath region I is a viscous layer, region II, which has an $O(R_{\unicode[STIX]{x1D6EC}}^{-1/2}\unicode[STIX]{x1D6EC})$ thickness. In this region, the centrifugal force is negligible, and so as in the flat-plate case, the streamwise and spanwise components of FSVD induce a quasi two-dimensional disturbance and a three-dimensional disturbance respectively, both having an $O(\unicode[STIX]{x1D716})$ amplitude (LWG). Like WZL, we focus on the three-dimensional disturbance, which will develop into Görtler vortices further downstream. In region III, which has an $O(\unicode[STIX]{x1D6EC})$ thickness, the perturbations are governed by steady or unsteady nonlinear boundary-region equations, but with the centrifugal force appearing in the wall-normal momentum equation. An outer region IV above region III, where turbulence undergoes equilibrium decay, has a transversal scale of $O(\unicode[STIX]{x1D6EC}R_{\unicode[STIX]{x1D6EC}})$ . With the assumption $r_{t}\ll 1$ , the disturbance remains linear in all the regions. When $r_{t}=O(1)$ , nonlinearity comes into play in regions III and IV (Wundrow & Goldstein Reference Wundrow and Goldstein2001). Görtler vortices then evolve nonlinearly, saturate and may undergo secondary instability.

In the boundary region (region III), where $\hat{x}$ and $y$ are both of $O(1)$ , the streamwise velocity has a magnitude greater than those of the normal and spanwise velocities by a factor of $O(R_{\unicode[STIX]{x1D6EC}})$ , while the pressure normalized by $\unicode[STIX]{x1D70C}U_{\infty }^{2}$ is of $O(R_{\unicode[STIX]{x1D6EC}}^{-2})$ . Therefore, we can write the velocity and pressure fields, $(u^{\ast },v^{\ast },w^{\ast })$ and $p^{\ast }$ , as

(2.8a,b ) $$\begin{eqnarray}(u^{\ast },v^{\ast },w^{\ast })/U_{\infty }=(u,R_{\unicode[STIX]{x1D6EC}}^{-1}v,R_{\unicode[STIX]{x1D6EC}}^{-1}w),\quad p^{\ast }/(\unicode[STIX]{x1D70C}U_{\infty }^{2})=R_{\unicode[STIX]{x1D6EC}}^{-2}p.\end{eqnarray}$$

Substitution of (2.4) and (2.8) with the Lamè coefficients, $(h_{1},h_{2},h_{3})=((r_{0}^{\ast }-y^{\ast })/r_{0}^{\ast },1,1)$ , into the Navier–Stokes (N–S) equations gives at leading-order equations for rescaled velocity $\boldsymbol{u}\equiv (u,v,w)$ and pressure $p$ (Hall Reference Hall1988),

(2.9a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{u}_{\hat{\unicode[STIX]{x1D70F}}}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}u^{2}\boldsymbol{j}=(0,-p_{y},p_{z})+(\unicode[STIX]{x2202}_{yy}^{2}+\unicode[STIX]{x2202}_{zz}^{2})\boldsymbol{u},\end{eqnarray}$$

where $\boldsymbol{j}$ is the unit vector in the wall-normal direction, the term containing $G_{\unicode[STIX]{x1D6EC}}$ reflects the essential influence of the wall curvature and $\unicode[STIX]{x1D712}^{-1}(\hat{x})$ is the scaled local radius of wall curvature. The system (2.9), which holds under the assumption of large curvature radius, will be referred to as the nonlinear unsteady boundary-region equations (NUBRE). A noteworthy property of the system is that it is parabolic in the streamwise direction but retains ellipticity in the spanwise direction.

The flow can be decomposed as a sum of the Blasius flow and the perturbation which arises as the response to FSVD, namely

(2.10) $$\begin{eqnarray}(u,v,w,p)=\left(U(x,\unicode[STIX]{x1D702}),V(x,\unicode[STIX]{x1D702}),0,-{\textstyle \frac{1}{2}}\right)+(\hat{u} ,\hat{v},{\hat{w}},\hat{p}).\end{eqnarray}$$

The Blasius boundary layer, taken as the unperturbed base flow, has the similarity solution,

(2.11a,b ) $$\begin{eqnarray}U=F^{\prime }(\unicode[STIX]{x1D702}),\quad V=(\unicode[STIX]{x1D702}F^{\prime }-F)/\sqrt{2\hat{x}}\quad \text{with}~\unicode[STIX]{x1D702}=y/\sqrt{2\hat{x}},\end{eqnarray}$$

where the prime denotes the differentiation with respect to $\unicode[STIX]{x1D702}$ , and the Blasius function $F$ is determined by the boundary-value problem,

(2.12a-d ) $$\begin{eqnarray}F^{\prime \prime \prime }+FF^{\prime \prime }=0,\quad F(0)=0,\quad F^{\prime }(0)=0,\quad F^{\prime }\rightarrow 1\quad \text{as}~\unicode[STIX]{x1D702}\rightarrow \infty .\end{eqnarray}$$

For $\unicode[STIX]{x1D702}\gg 1$ , the function $F\rightarrow \unicode[STIX]{x1D702}-\unicode[STIX]{x1D6FD}$ with $\unicode[STIX]{x1D6FD}\approx 1.217$ .

Inserting (2.10) into (2.9), we obtain the governing equations for the perturbation,

(2.13) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\hat{u} _{\hat{x}}+\hat{v}_{y}+{\hat{w}}_{z}=0,\\ \hat{u} _{\hat{\unicode[STIX]{x1D70F}}}+U\hat{u} _{\hat{x}}+\hat{u} U_{\hat{x}}+\hat{v}U_{y}+V\hat{u} _{y}=\hat{u} _{yy}+\hat{u} _{zz}+Q_{1},\\ \hat{v}_{\hat{\unicode[STIX]{x1D70F}}}+U\hat{v}_{\hat{x}}+\hat{u} V_{\hat{x}}+V\hat{v}_{y}+\hat{v}V_{y}+2G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}U\hat{u} =-\hat{p}_{y}+\hat{v}_{yy}+\hat{v}_{zz}+Q_{2},\\ {\hat{w}}_{\hat{\unicode[STIX]{x1D70F}}}+U{\hat{w}}_{\hat{x}}+V{\hat{w}}_{y}=-\hat{p}_{z}+{\hat{w}}_{yy}+{\hat{w}}_{zz}+Q_{3},\\ \end{array}\right\}\end{eqnarray}$$

where $Q_{1}$ , $Q_{2}$ and $Q_{3}$ represent the components of the nonlinear term $-(\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\hat{\boldsymbol{u}}-G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}\hat{u} ^{2}\boldsymbol{j}$ .

In the present nonlinear regime, the disturbance consists of all temporal and spanwise harmonics, and can be expressed as

(2.14) $$\begin{eqnarray}\displaystyle (\hat{u} ,\hat{v},{\hat{w}},\hat{p}) & = & \displaystyle r_{t}\mathop{\sum }_{m,n=-\infty }^{+\infty }((2\hat{x})\hat{u} _{m,n}(\hat{x},\unicode[STIX]{x1D702}),(2\hat{x})^{1/2}\hat{v}_{m,n}(\hat{x},\unicode[STIX]{x1D702}),{\hat{w}}_{m,n}(\hat{x},\unicode[STIX]{x1D702})/k_{3},\hat{p}_{m,n}(\hat{x},\unicode[STIX]{x1D702}))\nonumber\\ \displaystyle & & \displaystyle \times \,\exp (-\text{i}m\hat{k}_{1}\hat{\unicode[STIX]{x1D70F}}+\text{i}nk_{3}z).\end{eqnarray}$$

The factor $2\hat{x}$ in the streamwise velocity is introduced to offset the small divisor in numerical computations. The reality of the physical quality requires

(2.15) $$\begin{eqnarray}\hat{q}_{-m,-n}=(\hat{q}_{m,n})_{cc},\end{eqnarray}$$

where $\hat{q}_{m,n}$ stands for any of $(\hat{u} _{m,n},\hat{v}_{m,n},{\hat{w}}_{m,n},\hat{p}_{m,n})$ and the subscript cc indicates the complex conjugate. Substitution of (2.14) into (2.13) yields the equations for the Fourier coefficients as follows.

The continuity equation:

(2.16) $$\begin{eqnarray}\displaystyle 2\hat{u} _{m,n}+2\hat{x}\frac{\unicode[STIX]{x2202}\hat{u} _{m,n}}{\unicode[STIX]{x2202}\hat{x}}-\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}\hat{u} _{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x2202}\hat{v}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\text{i}n{\hat{w}}_{m,n}=0. & & \displaystyle\end{eqnarray}$$

The momentum equations:

(2.17) $$\begin{eqnarray}\displaystyle & & \displaystyle [2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})+2F^{\prime }-\unicode[STIX]{x1D702}F^{\prime \prime }]\hat{u} _{m,n}+2\hat{x}F^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}\hat{u} _{m,n}}{\unicode[STIX]{x2202}\hat{x}}}-F{\displaystyle \frac{\unicode[STIX]{x2202}\hat{u} _{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad -\,{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\hat{u} _{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}}+F^{\prime \prime }\hat{v}_{m,n}=-r_{t}(2\hat{x})\hat{l}_{m,n},\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & & \displaystyle [-\unicode[STIX]{x1D702}^{2}F^{\prime \prime }-\unicode[STIX]{x1D702}F^{\prime }+F+2(2\hat{x})^{3/2\;}G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}F^{\prime }]\hat{u} _{m,n}+[2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})+\!F^{\prime }+\unicode[STIX]{x1D702}F^{\prime \prime }]\hat{v}_{m,n}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,2\hat{x}F^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}\hat{v}_{m,n}}{\unicode[STIX]{x2202}\hat{x}}}-F{\displaystyle \frac{\unicode[STIX]{x2202}\hat{v}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}}-{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\hat{v}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\hat{p}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}}=-r_{t}(2\hat{x})\hat{e}_{m,n},\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & & \displaystyle [2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2}){\hat{w}}_{m,n}]+2\hat{x}F^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}{\hat{w}}_{m,n}}{\unicode[STIX]{x2202}\hat{x}}}-F^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}{\hat{w}}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad -\,{\displaystyle \frac{\unicode[STIX]{x2202}^{2}{\hat{w}}_{m,n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}}+2\hat{x}(\text{i}nk_{3}^{2})\hat{p}_{m,n}=-r_{t}(2\hat{x}){\hat{h}}_{m,n}.\end{eqnarray}$$

In the above equations, we have put

(2.20) $$\begin{eqnarray}\displaystyle \hat{l}_{m,n}=\left\{2\hat{x}\frac{\unicode[STIX]{x2202}(\widehat{u_{a}u_{a}})}{\unicode[STIX]{x2202}\hat{x}}+4\widehat{u_{a}u_{a}}-\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\widehat{u_{a}u_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x2202}(\widehat{u_{a}v_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{1}{k_{3}}\frac{\unicode[STIX]{x2202}(\widehat{u_{a}w_{a}})}{\unicode[STIX]{x2202}z}\right\}_{m,n}, & & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle \hat{e}_{m,n} & = & \displaystyle \left\{\left(2\hat{x}\right)\frac{\unicode[STIX]{x2202}(\widehat{v_{a}u_{a}})}{\unicode[STIX]{x2202}\hat{x}}+3\widehat{v_{a}u_{a}}-\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\widehat{v_{a}u_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{\unicode[STIX]{x2202}(\widehat{v_{a}v_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{1}{k_{3}}\frac{\unicode[STIX]{x2202}(\widehat{v_{a}w_{a}})}{\unicode[STIX]{x2202}z}+G_{\unicode[STIX]{x1D6EC}}\left(2\hat{x}\right)^{3/2}\unicode[STIX]{x1D712}\widehat{{u_{a}}^{2}}\right\}_{m,n},\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle {\hat{h}}_{m,n}=\left\{\left(2\hat{x}\right)\frac{\unicode[STIX]{x2202}(\widehat{w_{a}u_{a}})}{\unicode[STIX]{x2202}\hat{x}}+2\widehat{w_{a}u_{a}}-\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\widehat{u_{a}u_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x2202}(\widehat{w_{a}v_{a}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{1}{k_{3}}\frac{\unicode[STIX]{x2202}(\widehat{w_{a}w_{a}})}{\unicode[STIX]{x2202}z}\right\}_{m,n}. & & \displaystyle\end{eqnarray}$$

The new variables $(u_{a},v_{a},w_{a})$ are introduced for convenience of computation,

(2.23) $$\begin{eqnarray}(u_{a},v_{a},w_{a})=\mathop{\sum }_{m,n=-\infty }^{+\infty }(\hat{u} _{m,n}(\hat{x},\unicode[STIX]{x1D702}),\hat{v}_{m,n}(\hat{x},\unicode[STIX]{x1D702}),{\hat{w}}_{m,n}(\hat{x},\unicode[STIX]{x1D702}))\exp (-\text{i}m\hat{k}_{1}\hat{\unicode[STIX]{x1D70F}}+\text{i}nk_{3}z),\end{eqnarray}$$

and the $\widehat{\qquad }$ symbol in the equations denotes Fourier transform. For $r_{t}\ll 1$ , the nonlinear terms in the NUBRE are of $O(r_{t}^{2})$ , and so the LUBRE are recovered by neglecting those terms (LWG, WZL). For $r_{t}=O(1)$ , the nonlinear terms in the streamwise and spanwise momentum equations are exactly the same as those in the flat-plate nonlinear problem (Ricco et al. Reference Ricco, Luo and Wu2011), but the curved wall contributes an additional nonlinear terms $G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}\widehat{u_{a}^{2}}$ (i.e. $G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}\hat{u} ^{2}$ ) in the normal-wall momentum equation.

2.2 Initial and outer boundary conditions

Since the velocity fluctuation induced in the boundary layer is of small amplitude near the leading edge, nonlinearity is weak in this region, and hence the initial condition remains the same as that in the linear case. For the special case, where the oncoming disturbance consists of a pair of oblique components, the upstream conditions can be expressed as (LWG, WZL)

(2.24) $$\begin{eqnarray}\hat{u} _{1,\pm 1}\rightarrow q_{\pm }(U_{0}+(2\hat{x})^{1/2}U_{1}),\end{eqnarray}$$

(2.25) $$\begin{eqnarray}\displaystyle \hat{v}_{1,\pm 1} & \rightarrow & \displaystyle q_{\pm }\left\{V_{0}+(2\hat{x})^{1/2}V_{1}+\frac{\text{i}\text{e}^{\text{i}\hat{k}_{1}\hat{x}-\text{i}k_{2}(2\hat{x})^{1/2}\unicode[STIX]{x1D6FD}-(k_{3}^{2}+k_{2}^{2})\hat{x}}}{(k_{2}-\text{i}|k_{3}|)(2\hat{x})^{1/2}}[\text{e}^{\text{i}k_{2}(2\hat{x})^{1/2}\unicode[STIX]{x1D702}}-\text{e}^{-|k_{3}|(2\hat{x})^{1/2}\unicode[STIX]{x1D702}}]\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\text{e}^{-|k_{3}|(2\hat{x})^{1/2}\unicode[STIX]{x1D702}}\left(\frac{\unicode[STIX]{x1D6FD}}{4}+\sqrt{2\hat{x}}c_{1}\right)-\bar{v}_{c}\right\},\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle {\hat{w}}_{1,\pm 1} & \rightarrow & \displaystyle \mp \text{i}q_{\pm }\left\{W_{0}+(2\hat{x})^{1/2}W_{1}+{\displaystyle \frac{k_{2}\text{e}^{\text{i}\hat{k}_{1}\hat{x}-\text{i}k_{2}(2\hat{x})^{1/2}\unicode[STIX]{x1D6FD}-(k_{3}^{2}+k_{2}^{2})\hat{x}}}{(k_{2}-\text{i}|k_{3}|)}}[\text{e}^{\text{i}k_{2}(2\hat{x})^{1/2}\unicode[STIX]{x1D702}}-\text{e}^{-|k_{3}|(2\hat{x})^{1/2\unicode[STIX]{x1D702}}}]\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\text{e}^{-|k_{3}|(2\hat{x})^{1/2}\unicode[STIX]{x1D702}}\left[1-\unicode[STIX]{x1D6FD}\left(\text{i}(k_{2}+\text{i}|k_{3}|)+\frac{3}{4}|k_{3}|\right)\sqrt{2\hat{x}}\right]-\bar{w}_{c}\right\},\end{eqnarray}$$

as $\hat{x}\rightarrow 0$ , where $q_{\pm }=\pm \text{i}k_{3}(\hat{u} _{3,\pm }^{\infty }\pm \text{i}\hat{u} _{2,\pm }^{\infty })$ ,

(2.27) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{v}_{c}=-\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D6FD}}{4}+\left(2\hat{x}\right)^{1/2}\left\{-\frac{\text{i}}{2}\left(k_{2}+\text{i}\left|k_{3}\right|\right)\unicode[STIX]{x1D702}^{2}+\unicode[STIX]{x1D6FD}\left(\text{i}k_{2}-\frac{1}{4}|k_{3}|\right)\unicode[STIX]{x1D702}+c_{1}\right\}, & \displaystyle\end{eqnarray}$$

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{w}_{c}=1+\left(2\hat{x}\right)^{1/2}\left\{\text{i}\left(k_{2}+\text{i}|k_{3}|\right)\unicode[STIX]{x1D702}-\unicode[STIX]{x1D6FD}\left(\text{i}k_{2}-{\textstyle \frac{1}{4}}|k_{3}|\right)\right\}. & \displaystyle\end{eqnarray}$$

The constant $c_{1}$ is obtained numerically along with $U_{k},V_{k}$ and $W_{k}$ ( $k=0,1$ ) by solving the system (B 1)–(B 8) on page 200 in WZL. As in LWG and WZL, here we neglected the so-called ‘quasi two-dimensional’ particular solution as well as the homogeneous solution whose streamwise velocity behaves like $x^{\unicode[STIX]{x1D70E}}$ ( $\unicode[STIX]{x1D70E}\approx 0.213$ ). The latter is not excited by FSVD in the present context, but its nonlinear development (with an arbitrarily prescribed amplitude) and possible role in transition was studied by Martin et al. (Reference Martin, Martel, Paredes and Theofilis2015) and references therein.

The far-field condition is derived by considering the outer region IV. Similar to the flat-plate cases investigated by Wundrow & Goldstein (Reference Wundrow and Goldstein2001) and Ricco et al. (Reference Ricco, Luo and Wu2011), the perturbation in region IV is influenced by the viscous motion in region III through the displacement effect. It was shown there that the normal velocity has the asymptote

(2.29) $$\begin{eqnarray}v\rightarrow \bar{\unicode[STIX]{x1D6FF}}_{\hat{x}}(\hat{x},\hat{\unicode[STIX]{x1D70F}})\quad \text{as}~y\rightarrow \infty ,\end{eqnarray}$$

where $\bar{\unicode[STIX]{x1D6FF}}$ is the spanwise-averaged displacement thickness,

(2.30) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6FF}}(\hat{x},\hat{\unicode[STIX]{x1D70F}})=\frac{k_{3}}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}/k_{3}}\int _{0}^{\infty }[1-u(y,z,\hat{x},\hat{\unicode[STIX]{x1D70F}})]\,\text{d}y\,\text{d}z.\end{eqnarray}$$

The analysis of the outer region IV and its interaction with region III is presented in appendix A. Matching with the outer solution gives the far-field condition,

(2.31) $$\begin{eqnarray}(\hat{u} _{m,n},\hat{v}_{m,n},{\hat{w}}_{m,n},\hat{p}_{m,n})\rightarrow (0,(2\hat{x})^{-1/2}v_{m,n}^{\dagger },k_{3}w_{m,n}^{\dagger },\unicode[STIX]{x1D716}R_{\unicode[STIX]{x1D6EC}}p_{m,n}^{\dagger })\quad \text{as}~\unicode[STIX]{x1D702}\rightarrow \infty ,\end{eqnarray}$$

for $n\neq 0$ , where $v_{m,n}^{\dagger },w_{m,n}^{\dagger }$ and $p_{m,n}^{\dagger }$ are given in appendix A.

For steady perturbations, the initial and far-field conditions can be further simplified. Since $\hat{k}_{1}$ is zero, only spanwise spectral components exist, and so $(\hat{u} _{-n},\hat{v}_{-n},{\hat{w}}_{-n},\hat{p}_{-n})=(\hat{u} _{n},\hat{u} _{n},{\hat{w}}_{n},\hat{p}_{n})_{cc}$ . The boundary conditions are matched with the outer solution in the same way as in the unsteady case.

The initial-boundary-value problem, consisting of (2.17)–(2.20), (2.25)–(2.27) and (2.31), governs the excitation, linear and nonlinear evolution of Görtler vortices in the presence of FSVD for the generic case of $G_{\unicode[STIX]{x1D6EC}}=O(1)$ . Several features of the formulation are worth emphasizing. First, the initial (upstream) and boundary conditions are interrelated, both being determined by FSVD. There is inherent consistency between them: the upstream limit ( $\hat{x}\rightarrow 0$ ) of the latter matches the large- $y$ limit of the former. Second, throughout the entire course of development, including the fully developed phase when their streamwise velocity has acquired an $O(1)$ magnitude, Görtler vortices are continuingly influenced by FSVD via the forcing of the transverse velocities, that is, the boundary condition is as an inseparable part of the evolution just as the initial condition is. Finally, the viscous motion in the boundary layer is not a mere passive response to FSVD; it affects simultaneously the disturbances in the outer region as well.

3.1 Discretized equations

The initial-boundary-value problem, equations (2.17)–(2.20) supplemented by (2.25)–(2.27) and (2.31), is to be solved numerically by using appropriate finite-difference schemes.

For the purpose of computation, we define $\hat{f}=\hat{u} _{\unicode[STIX]{x1D702}}$ and ${\hat{g}}={\hat{w}}_{\unicode[STIX]{x1D702}}$ , and recast equations (2.17)–(2.20) into a system of first-order equations,

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\hat{u} _{\unicode[STIX]{x1D702}}=\hat{f},\quad \hat{v}_{\unicode[STIX]{x1D702}}=-2\hat{u} -\text{i}n{\hat{w}}+\unicode[STIX]{x1D702}\hat{f}-2\hat{x}\hat{u} _{\hat{x}},\quad {\hat{w}}_{\unicode[STIX]{x1D702}}={\hat{g}},\\ \hat{f}_{\unicode[STIX]{x1D702}}=[2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})+2F^{\prime }-\unicode[STIX]{x1D702}F^{\prime \prime }]\hat{u} +F^{\prime \prime }\hat{v}-F\hat{f}+2\hat{x}F^{\prime }\hat{u} _{\hat{x}}-\hat{l}_{m,n},\\ \hat{p}_{\unicode[STIX]{x1D702}}=[3\unicode[STIX]{x1D702}F^{\prime }-3F-2(2\hat{x})^{3/2}G_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D712}F^{\prime }+\unicode[STIX]{x1D702}2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})]\hat{u} \\ +\,[-2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})-F^{\prime }]\hat{v}-\text{i}nF{\hat{w}}-\hat{f}-\text{i}n{\hat{g}}\\ +\,2\hat{x}(-F+\unicode[STIX]{x1D702}F^{\prime })\hat{u} _{\hat{x}}-2\hat{x}F^{\prime }\hat{v}_{\hat{x}}-2\hat{x}\hat{f}_{\hat{x}}-\unicode[STIX]{x1D702}\hat{l}_{m,n}+\hat{e}_{m,n},\\ {\hat{g}}_{\unicode[STIX]{x1D702}}=[2\hat{x}(-\text{i}m\hat{k}_{1}+n^{2}k_{3}^{2})]{\hat{w}}+2\hat{x}(\text{i}nk_{3}^{2})\hat{p}-F^{\prime }{\hat{g}}+2\hat{x}F^{\prime }{\hat{w}}_{\hat{x}}-{\hat{h}}_{m,n},\end{array}\right\}\end{eqnarray}$$

which can be rewritten in a vector form,

(3.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=C_{0}\boldsymbol{U}+C_{1}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}\hat{x}}+f(\boldsymbol{U}),\end{eqnarray}$$

where $\boldsymbol{U}=(\hat{u} ,\hat{v},{\hat{w}},\hat{f},\hat{p},{\hat{g}})$ is the solution vector, $C_{0}$ and $C_{1}$ are the coefficient matrices and $f(\boldsymbol{U})$ stands for nonlinear terms. For brevity, the subscripts ‘ $m,n$ ’ are omitted.

For components with $n=0$ , $\hat{u}$ and $\hat{v}$ are decoupled from ${\hat{w}}$ , and the system (3.1) can be simplified to three equations for $\hat{u} ,\hat{v}$ and $\hat{f}$ :

(3.3) $$\begin{eqnarray}\left.\begin{array}{@{}cc@{}} & \hat{u} _{\unicode[STIX]{x1D702}}=\hat{f},\quad \hat{v}_{\unicode[STIX]{x1D702}}=-2\hat{u} +\unicode[STIX]{x1D702}\hat{f}-2\hat{x}\hat{u} _{\hat{x}},\\ & \hat{f}_{\unicode[STIX]{x1D702}}=[2\hat{x}(-\text{i}m\hat{k}_{1})+2F^{\prime }-\unicode[STIX]{x1D702}F^{\prime \prime }]\hat{u} +F^{\prime \prime }\hat{v}-F\hat{f}+2\hat{x}F^{\prime }\hat{u} _{\hat{x}}-\hat{l}_{m,0};\end{array}\right\}\end{eqnarray}$$

and two equations for ${\hat{w}}$ and ${\hat{g}}$ ,

(3.4a,b ) $$\begin{eqnarray}{\hat{w}}_{\unicode[STIX]{x1D702}}={\hat{g}},\quad {\hat{g}}_{\unicode[STIX]{x1D702}}=2\hat{x}(-\text{i}m\hat{k}_{1}){\hat{w}}-F^{\prime }{\hat{g}}+2\hat{x}F^{\prime }{\hat{w}}_{\hat{x}}-{\hat{h}}_{m,0}.\end{eqnarray}$$

The NUBRE are parabolic in the streamwise direction, and hence can be solved by a marching procedure in $\hat{x}$ -direction. The Malik scheme (Malik Reference Malik1990) is employed in $\unicode[STIX]{x1D702}$ -direction, and the second-order backward finite-difference scheme is applied in the streamwise direction, namely

(3.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{U}_{j+1}-\boldsymbol{U}_{j}=\frac{h}{2}\left({\boldsymbol{U}^{\prime }}_{j+1}+{\boldsymbol{U}^{\prime }}_{j}\right),\\ \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{U}^{n+1}}{\unicode[STIX]{x2202}\hat{x}}=\left.\left[\frac{3}{2}\boldsymbol{U}^{n+1}-2\boldsymbol{U}^{n}+\frac{1}{2}\boldsymbol{U}^{n-1}\right]\right/\unicode[STIX]{x0394}\hat{x},\end{array}\right\}\end{eqnarray}$$

where $j$ and $n$ indicate the position of mesh points in the transverse and streamwise directions respectively. Grouping all terms at the streamwise location $n+1$ on the left-hand side (except the nonlinear terms), we end up with the discrete equations,

(3.6) $$\begin{eqnarray}\left[\unicode[STIX]{x1D644}-\frac{h}{2}c_{j+1}\right]\boldsymbol{U}_{j+1}^{n+1}+\left[-\unicode[STIX]{x1D644}-\frac{h}{2}c_{j}\right]\boldsymbol{U}_{j}^{n+1}=\frac{h}{2}\left[g_{j+1}+g_{j}\right]+\frac{h}{2}[f(\boldsymbol{U}_{j+1}^{n+1})+f(\boldsymbol{U}_{j}^{n+1})],\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the unit matrix, and we have put

(3.7) $$\begin{eqnarray}\displaystyle f(\boldsymbol{U}_{j}^{n+1})=\left(\begin{array}{@{}cccccc@{}}0, & 0, & 0, & -\hat{l}, & -\unicode[STIX]{x1D702}\hat{l}+\hat{e}, & -{\hat{h}}\end{array}\right)_{j}^{n+1}, & & \displaystyle\end{eqnarray}$$

(3.8a,b ) $$\begin{eqnarray}\displaystyle c_{j+1}=C_{_{0,j+1}}+{\displaystyle \frac{3}{2\unicode[STIX]{x0394}\hat{x}}}C_{1,j+1},\quad g_{j+1}=C_{1,j+1}\left.\left(-2\boldsymbol{U}_{j+1}^{n}+\frac{1}{2}\boldsymbol{U}_{j+1}^{n-1}\right)\right/\unicode[STIX]{x0394}\hat{x}. & & \displaystyle\end{eqnarray}$$

On introducing the discrete solution vector,

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{n+1}=\left(\begin{array}{@{}ccccc@{}}\boldsymbol{U}_{J}^{n+1}, & \boldsymbol{U}_{J-1}^{n+1}, & \ldots , & \boldsymbol{U}_{1}^{n+1}, & \boldsymbol{U}_{0}^{n+1}\end{array}\right)^{\text{T}},\end{eqnarray}$$

with $J$ being the number of mesh points in the wall-normal direction, the system of nonlinear algebraic equations, (3.6), is written in the matrix form,

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D647}\unicode[STIX]{x1D6F7}_{n+1}=\boldsymbol{A}+\boldsymbol{N}_{n+1},\end{eqnarray}$$

where $\unicode[STIX]{x1D647}$ is the linear operator, $\boldsymbol{A}$ stands for the terms whose values are evaluated at upstream positions and $\boldsymbol{N}_{n+1}$ denotes the nonlinear terms (see appendix B for their expressions). Only half of the Fourier modes needs to be computed since Fourier modes with negative indices $m$ can be evaluated by (2.15). The nonlinear terms are evaluated by applying the pseudo-spectral method, that is, the velocities in spectral space are transformed back to physical space, where multiplications are carried out, and then those products are Fourier transformed to obtain the nonlinear terms in spectral space. In order to prevent aliasing error, i.e. the spurious energy cascade from unresolved high-frequency modes into the resolved low-frequency ones, the 3/2-rule is applied (Kim, Moin & Moser Reference Kim, Moin and Moser1987). Seventeen Fourier modes are retained for capturing nonlinear effects. The domain size in the $\unicode[STIX]{x1D702}$ -direction is 30, within which 1000 grid points are deployed.

3.2 Solution procedure and validation of the algorithm

With the nonlinear terms being treated implicitly, the system (3.10) consists of nonlinear algebraic equations. It will be solved by using an iterative procedure involving a predictor–corrector and underrelaxation.

First, the nonlinear terms are evaluated by using the value at the previous streamwise point. This leads to a linear system of simultaneous equations,

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D647}\unicode[STIX]{x1D6F7}_{n+1}^{\ast }=\boldsymbol{A}+\boldsymbol{N}_{n},\end{eqnarray}$$

in which the coefficient matrix $\unicode[STIX]{x1D647}$ is block tri-diagonal. This explicit method works when $G_{\unicode[STIX]{x1D6EC}}<20$ according to our computation, but fails to converge and a non-physical result appears when $G_{\unicode[STIX]{x1D6EC}}>20$ .

An iteration method is thus applied in order to obtain the convergent solution. The $\unicode[STIX]{x1D6F7}_{n+1}^{\ast }$ predicted by the explicit scheme is used to re-evaluate the nonlinear terms. Then we solve the resulting equations to get a new value of $\unicode[STIX]{x1D6F7}_{n+1}$ , which we use to compute $\boldsymbol{N}_{n+1}$ . This step is repeated, and the correct nonlinear terms and the solution are obtained when the difference between the two consecutive outcomes for $\unicode[STIX]{x1D6F7}_{n+1}$ is less than a prescribed tolerance at all grid points. The non-physical result does not appear any more in the result, which is shown in symbols in figure 2(a). However, this iterative method is feasible only up to a certain streamwise position, beyond which it failed to converge, indicating that the predicted nonlinear terms are not suitable when the amplitude is sufficiently large and nonlinearity is strong.

Figure 2. Comparison of the downstream development of the maximum streamwise velocity of different Fourier modes for $\unicode[STIX]{x1D716}=0.0001$ . (a) The results computed without the underrelaxation iteration (symbols) and with the underrelaxtion iteration (solid lines) for $G_{\unicode[STIX]{x1D6EC}}=1344$ and $R_{\unicode[STIX]{x1D6EC}}=1222$ . (b) The results computed by using streamwise step $\unicode[STIX]{x0394}\hat{x}=0.005$ (symbols) and $\unicode[STIX]{x0394}\hat{x}=0.001$ (solid lines) for $G_{\unicode[STIX]{x1D6EC}}=1765$ and $R_{\unicode[STIX]{x1D6EC}}=2298$ . The vertical dashed lines indicate the last streamwise position before the calculation without underrelaxation blows up.

Frequently, the convergence to a solution can be enhanced by a technique called underrelaxation, which was also used in solving the so-called parabolized stability equations (Zhao et al. Reference Zhao, Zhang, Liu and Luo2016). We use the value of $\boldsymbol{N}_{n}$ obtained at the previous location, and $\boldsymbol{N}_{n+1}^{\ast }$ obtained by using $\unicode[STIX]{x1D6F7}_{n+1}^{\ast }$ , to extrapolate a value for $\hat{\boldsymbol{N}}_{n+1}$ as follows:

(3.12) $$\begin{eqnarray}\hat{\boldsymbol{N}}_{n+1}=\boldsymbol{N}_{n}+\unicode[STIX]{x1D71B}(\boldsymbol{N}_{n+1}^{\ast }-\boldsymbol{N}_{n}),\end{eqnarray}$$

where $0<\unicode[STIX]{x1D71B}<1$ is the parameter controlling underrelaxation, and $\hat{\boldsymbol{N}}_{n+1}$ is the new nonlinear terms for the next iterative step. Hence, our predictor–corrector procedure follows the steps below.

Step 1: The nonlinear terms $\boldsymbol{N}_{n}$ is evaluated by using the value of $\unicode[STIX]{x1D731}_{n}$ at the last streamwise location.

Step 2: An estimated value of $\unicode[STIX]{x1D731}_{n+1}^{(i)}$ is obtained by solving the tri-diagonal form of the equations, and $\unicode[STIX]{x1D731}_{n+1}^{(i)}$ is utilized to obtain upgraded nonlinear terms $\boldsymbol{N}_{n+1}^{(i)}$ , where $i$ refers to the number of iterations. Using $\boldsymbol{N}_{(n+1)}^{(i+1)}=\boldsymbol{N}_{n}+\unicode[STIX]{x1D71B}(\boldsymbol{N}_{n+1}^{(i)}-\boldsymbol{N}_{n})$ , we can update the value of $\unicode[STIX]{x1D731}_{n+1}$ , which is written as $\unicode[STIX]{x1D731}_{n+1}^{(i+1)}$ .

Step 3: Repeat step 2 until the amplitude of every Fourier mode converges. A satisfactory convergence is achieved rather rapidly by this method by choosing $\unicode[STIX]{x1D71B}=0.2$ .

Different streamwise step lengths are also chosen to ensure convergence and accuracy of the solutions. Our numerical results indicate that the streamwise step needs to be reduced when the Görtler number is large. For instance, when $\unicode[STIX]{x0394}\hat{x}=0.005$ the computation converges for $G_{\unicode[STIX]{x1D6EC}}=1344$ as shown in figure 2(a), but diverges for $G_{\unicode[STIX]{x1D6EC}}=1765$ , which requires a smaller step size $\unicode[STIX]{x0394}\hat{x}=0.001$ . A comparison of the results using these two step sizes for $G_{\unicode[STIX]{x1D6EC}}=1765$ is displayed in figure 2(b). In all our calculations, we set $\unicode[STIX]{x1D712}=1$ since the curved plates in experiments often have a constant curvature.

Figure 3. Validation of numerical algorithm and code. (a) Comparison of the downstream development of the maximum streamwise velocity of Görtler vortices with the linear results of WZL. The parameters are $G_{\unicode[STIX]{x1D6EC}}=1804$ , $R_{\unicode[STIX]{x1D6EC}}=2192,\unicode[STIX]{x1D716}=7.0\times 10^{-7}$ and $G_{\unicode[STIX]{x1D6EC}}=1344,R_{\unicode[STIX]{x1D6EC}}=1222,\unicode[STIX]{x1D716}=7.0\times 10^{-6}$ . (b) Comparison of the nonlinear development of $\max u_{rms}^{\prime }$ , the maximum root mean square of the time-dependent streamwise velocity of streaks (symbols), with the results of Ricco et al. (Reference Ricco, Luo and Wu2011) (lines) for fixed $\hat{k}_{1}=2$ , but different $\unicode[STIX]{x1D716}=0.005$ (triangle and solid line), 0.01 (circle and dashed lines) and 0.015 (square and dash-dotted lines).

In order to validate our code, two tests are carried out.

1. (i) Our code is designed primarily for nonlinear problems, but it must be able to reproduce the linear solution when $r_{t}\ll 1$ . The results obtained by the present method are compared in figure 3(a) with the linear results of WZL, shown in figure 11 of their paper. The parameters adopted are $R_{\unicode[STIX]{x1D6EC}}=1804,G_{\unicode[STIX]{x1D6EC}}=2192$ and $\unicode[STIX]{x1D716}=7.0\times 10^{-7}$ , and $R_{\unicode[STIX]{x1D6EC}}=1344$ , $G_{\unicode[STIX]{x1D6EC}}=1222$ and $\unicode[STIX]{x1D716}=7.0\times 10^{-6}$ . The linear result was rescaled to fit the nonlinear result at the first location. Excellent agreement is seen for the entire development of the maximum streamwise velocity as well as for the growth rate of steady vortices.

2. (ii) Our code pertains to a boundary layer over a concave wall with $G_{\unicode[STIX]{x1D6EC}}>0$ . When $G_{\unicode[STIX]{x1D6EC}}=0$ , it should be able to reproduce the results for the flat-plate case, which was computed by Ricco et al. (Reference Ricco, Luo and Wu2011) using a scheme entirely different from ours. The intensity of streaks or Görtler vortices is represented by the root mean square (r.m.s.) of the streamwise velocity, defined as (see Pope 2000, p. 687)

   (3.13) $$\begin{eqnarray}u_{rms}\equiv r_{t}\left[\mathop{\sum }_{m,n}|2\hat{x}\hat{u} _{m,n}|^{2}\right]^{1/2}.\end{eqnarray}$$
   In order to measure time-dependent components, we also introduce $u_{rms}^{\prime }$ , which has the same expression as $u_{rms}$ except that the zero-frequency components are excluded in the summation. The downstream development of $\max u_{rms}^{\prime }$ for three different values of $\unicode[STIX]{x1D716}$ is compared in figure 3(b), and an excellent agreement is seen for all three cases. This gives us the confidence that our code can correctly predict the nonlinear evolution of streaks or Görtler vortices.

4.1 Comparison with Swearingen and Blackweld (1987)

SB studied the growth, breakdown and transition to turbulence of Görtler vortices by means of smoke-wire visualization and multiple-probe hot-wire measurements in a low-speed wind tunnel. The free-stream velocity $U_{\infty }=5~\text{ms}^{-1}$ and the boundary layer developed over a concave wall with a radius of curvature $r_{0}^{\ast }=3.20$ m. The spanwise wavelength of the Görtler vortices was $\unicode[STIX]{x1D706}_{z}\approx 2.3$  cm. These values are used to calculate the non-dimensional parameters, $R_{\unicode[STIX]{x1D6EC}}=U_{\infty }\unicode[STIX]{x1D706}_{z}/(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D708})=1145$ and $G_{\unicode[STIX]{x1D6EC}}=R_{\unicode[STIX]{x1D6EC}}^{2}\unicode[STIX]{x1D706}_{z}/(2\unicode[STIX]{x03C0}r_{0}^{\ast })=1501$ , where the kinematic viscosity coefficient of dry air $\unicode[STIX]{x1D708}=1.6\times 10^{-5}~\text{m}^{2}~\text{s}^{-1}$ is taken. The non-dimensional streamwise distance $\hat{x}$ is calculated by the relation

(4.1) $$\begin{eqnarray}\hat{x}=x^{\ast }/\unicode[STIX]{x1D6EC}R_{\unicode[STIX]{x1D6EC}}=2\unicode[STIX]{x03C0}x^{\ast }/(R_{\unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D706}_{z}).\end{eqnarray}$$

The corresponding Reynolds and Görtler numbers as well as the non-dimensional streamwise locations of measurement are listed in table 1. We focus on the steady components $k_{1}=0$ , and take $\hat{u} _{3,\pm }^{\infty }=\pm 1$ and $k_{2}=k_{3}=1$ , since the perturbation is likely to be isotropic. It follows that $\hat{u} _{2,\pm 1}^{\infty }=-1.0$ . The computation domain starts from $\hat{x}=0.001$ .

Figure 4. Comparison of the downstream development of $\max u_{rms}$ for (a) fixed $G_{\unicode[STIX]{x1D6EC}}$ but different values of $\unicode[STIX]{x1D716}$ and (b) fixed $\unicode[STIX]{x1D716}$ but different values of $G_{\unicode[STIX]{x1D6EC}}$ .

The two important parameters are $\unicode[STIX]{x1D716}$ and $G_{\unicode[STIX]{x1D6EC}}$ . In order to demonstrate their role, we solved the initial-boundary-value problems for several $\unicode[STIX]{x1D716}$ and $G_{\unicode[STIX]{x1D6EC}}$ . We choose to measure the intensity of vortices by $u_{rms}$ , the r.m.s of all harmonic components (as defined by (3.13) but with the double summation reducing to a single sum for the present steady limit). Figure 4(a) shows the downstream development of $\max u_{rms}$ for $G_{\unicode[STIX]{x1D6EC}}=1501$ and four different values of $\unicode[STIX]{x1D716}$ . The flat-plate case ( $G_{\unicode[STIX]{x1D6EC}}=0$ ) is also computed for comparison. When $\hat{x}\rightarrow 0$ , the centrifugal force vanishes, and so the solutions for $G_{\unicode[STIX]{x1D6EC}}=1501$ and $G_{\unicode[STIX]{x1D6EC}}=0$ overlap for each value of $\unicode[STIX]{x1D716}$ . The centrifugal effect influences the evolution of the disturbance weakly in the range of $\hat{x}=0.03-0.05$ . As $\unicode[STIX]{x1D716}$ increases, the response to the FSVD becomes stronger, and Görtler vortices and streaks saturate more quickly. For the case of strongest FSVD, $\unicode[STIX]{x1D716}=0.035$ (which is characteristic of turbomachinery flow conditions), the linear amplification stage of Görtler vortices almost disappears, which means that it would be completely inappropriate to use local eigenmodes as initial conditions for integrating NUBRE. Interestingly, in the case of $G_{\unicode[STIX]{x1D6EC}}\neq 0$ , $u_{rms}$ eventually saturates at almost the same level for different values of $\unicode[STIX]{x1D716}$ , suggesting that the saturated state is an intrinsic property of the boundary-layer flow. This is however not the case for the flat-plate boundary layer. The responses of boundary layer to FSVD with $\unicode[STIX]{x1D716}=0.0007$ are compared for different $G_{\unicode[STIX]{x1D6EC}}$ in figure 4(b). It is found that $u_{rms}$ with a larger $G_{\unicode[STIX]{x1D6EC}}$ has a greater linear growth rate as expected, and saturates at a higher level. In the limit of large Görtler number, an asymptotic theory for the fully nonlinear stage was proposed by Hall & Lakin (Reference Hall and Lakin1988), but a direct comparison is beyond the scope of the present study.

Figure 5. The downstream development of $\max u_{1,rms}$ , the r.m.s. of the spanwise-dependent streamwise velocity: the solid line represents the theoretical prediction; the squares symbols indicate the experimental measurement of the steady components (SB), whereas the circle symbols indicate the unsteady contribution.

The intensity of FSVD is controlled by $\unicode[STIX]{x1D716}$ , but an alternative and more common measure is turbulence level $Tu$ . For the assumed form of perturbations, equation (2.5), the two are related by the equation: $Tu=2\unicode[STIX]{x1D716}(\hat{u} _{1}^{\infty 2}+\hat{u} _{2}^{\infty 2}+\hat{u} _{3}^{\infty 2})^{1/2}/\sqrt{3}$ . Our theory shows that the streamwise velocity $\hat{u} _{1}^{\infty }$ does not play a leading-order role, suggesting that a modified form involving only the transverse velocities, $Tu^{\dagger }=2\unicode[STIX]{x1D716}(\hat{u} _{2}^{\infty 2}+\hat{u} _{3}^{\infty 2})^{1/2}/\sqrt{2}$ , is a more appropriate characterization. We choose $\unicode[STIX]{x1D716}=0.0007$ in order to fit the measured streamwise velocity at the first location of measurement, and correspondingly $Tu^{\dagger }=0.14\,\%$ , and $Tu=0.14\,\%$ as well if one takes $\hat{u} _{1}^{\infty }=1$ . In the experiment of SB, the FTS level was lower than $0.07\,\%$ , but it most likely refers to unsteady fluctuations rather than steady components. Indeed SB stated that transverse variations of the steady flow were less than $0.5\,\%$ . We believe that it is this part of the perturbation that excites Görtler vortices, but further detailed experiments are needed to characterize it.

Figure 6. Comparison of iso- $U$ contours at different streamwise locations: (a) theoretical prediction; (b) measurements of SB. The iso- $U$ contours are plotted with an increment of 0.10 $U_{\infty }$ as were the measurements of SB.

As in the experiment of SB, we monitor the r.m.s. of the spanwise-dependent components in the streamwise velocity, defined as

(4.2) $$\begin{eqnarray}u_{1,rms}\equiv r_{t}\left[\mathop{\sum }_{n\neq 0}|2\hat{x}\hat{u} _{n}|^{2}\right]^{1/2}.\end{eqnarray}$$

Note that the spanwise uniform mean-flow distortion is excluded. We compare the development of the computed $u_{1,rms}$ with the measurements of SB in figure 5. A perfect agreement is obtained up to $\hat{x}=0.24$ . The r.m.s. of the fluctuations experiences an initial non-modal growth and then amplifies exponentially as indicated by the logarithmic scale of the figure. Nonlinear interactions take effect from approximately $\hat{x}\approx 0.17$ . Disagreement is noted after $\hat{x}=0.26$ , where secondary instability of Görtler vortices may have sufficiently developed, and high-frequency fluctuations become strong enough to interact nonlinearly to affect the evolution of the primary vortices (Girgis & Liu Reference Girgis and Liu2006). Disturbances in the experiments of SB contained time-dependent fluctuations. The unsteady contribution was obtained by subtracting the steady signals from the measured spanwise mean. The r.m.s. of unsteady fluctuations increases with $\hat{x}$ , but remains much smaller than that of steady vortices before $\hat{x}=0.215$ ( $x^{\ast }=90$  cm). Unsteady perturbations apparently experienced fast exponential growth in the range of $\hat{x}=0.22$ – $0.26$ , as indicated by the dashed line in figure 5. This issue will be discussed in connection with secondary instability of Görtler vortices.

Of particular interest is the streamwise velocity of the distorted boundary layer,

(4.3) $$\begin{eqnarray}U(\unicode[STIX]{x1D702},z;\hat{x},\hat{\unicode[STIX]{x1D70F}})\equiv F^{\prime }(\unicode[STIX]{x1D702})+r_{r}(2\hat{x})\mathop{\sum }_{m,n}\hat{u} _{m,n}(\hat{x},\unicode[STIX]{x1D702})\exp (-\text{i}m\hat{k}_{1}\hat{\unicode[STIX]{x1D70F}}+\text{i}nk_{3}z).\end{eqnarray}$$

Contours of the predicted and measured $U$ in the cross-section (the $y$ – $z$ plane) are displayed in figure 6 for comparison. The similarity variable $\unicode[STIX]{x1D702}$ is related to the dimensional coordinate $y^{\ast }$ via the relation

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D702}=y^{\ast }(U_{\infty }/2\unicode[STIX]{x1D708}x^{\ast })^{1/2}.\end{eqnarray}$$

The contours evolved into a mushroom-like structure, as was predicted by previous calculations (Hall Reference Hall1988; Lee & Liu Reference Lee and Liu1992; Li & Malik Reference Li and Malik1995). The present comparisons show that not only the intensity but also the shape of Görtler vortices are in good agreement up to $\hat{x}=0.238$ ( $x^{\ast }=100$  cm). The gross feature of vortices remains similar for $\hat{x}\geqslant 0.262$ ( $x^{\ast }\geqslant 110$  cm), but the measured contours appear smoother. This is probably an indication that secondary instability has truly developed by $\hat{x}=0.262$ with the resulting high-frequency fluctuations affecting Görtler vortices through the Reynolds stresses (Girgis & Liu Reference Girgis and Liu2006).

Figure 7 displays the profiles of the streamwise velocity at different downstream locations and the two spanwise positions corresponding to a peak and a valley. A valley/peak refers to a low-/high-speed region, where the velocity exhibits a deficit/excess from the Blasius flow, respectively, We make comparison at these positions since only limited data were given by SB. There is a perfect agreement at every streamwise measurement location except at the valley of $\hat{x}=0.191$ ( $x^{\ast }=80$  cm). The reason for this discrepancy remains unclear.

Figure 7. Comparison of the predicted streamwise velocity profiles with the experimental measurements (SB 1987) at the spanwise positions corresponding to a valley (a) and a peak (b). Solid lines, theoretical predictions; symbols, experimental results. The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501$ and $R_{\unicode[STIX]{x1D6EC}}=1145$ .

Figure 8. Comparisons of the predicted displacement thickness (a) and the skin friction (b) with the experimental measurements (SB 1987). Lines: theoretical results; circles: experimental measurement in the low-speed region; squares: experimental measurements in the high-speed region. The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501$ and $R_{\unicode[STIX]{x1D6EC}}=1145$ .

Due to the presence of Görtler vortices, the local boundary-layer thickness differs from that of the Blasius flow and becomes spanwise dependent. The local displacement thickness is defined as

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{\ast }(x,z)=\int _{0}^{\infty }[1-u(x,y,z)]\,\text{d}y^{\ast }=\left(\frac{2\unicode[STIX]{x1D708}x^{\ast }}{U_{\infty }}\right)^{1/2}\int _{0}^{\infty }[1-u(x,\unicode[STIX]{x1D702},z)]\,\text{d}\unicode[STIX]{x1D702}.\end{eqnarray}$$

Figure 8(a) compares the predicted displacement thickness $\unicode[STIX]{x1D6FF}^{\ast }$ with the experimental measurements of SB in the high-speed (peak) and low-speed (valley) regions. From $x^{\ast }=20$  cm to 80 cm roughly, the displacement thickness in the low-speed region increases quickly while that in the high-speed region decreases, both because of exponential amplification of the vortices. The variations (increase and decrease) become moderate from $x^{\ast }=80$  cm to $90$  cm, due to saturation of Görtler vortices. The displacement thicknesses in the high- and low-speed regions are both in good agreement with the experimental data from $x^{\ast }=10$  cm to $90$  cm. The difference between the prediction and the measurement in the low-speed region near $x^{\ast }=95$  cm may be attributed to inaccuracy of measurements since this difference almost disappears in the high-speed region. The displacement thickness in the low-speed region decreases slightly at $x^{\ast }=100$  cm because of the over saturation of Görtler vortices, as indicated by figure 5. The comparison diverges downstream of $x^{\ast }=105$  cm, where secondary instability takes over and influences the displacement thickness substantially. Figure 8(b) displays the predicted skin friction values $\unicode[STIX]{x1D70F}^{\ast }=(\unicode[STIX]{x2202}u^{\ast }/\unicode[STIX]{x2202}y^{\ast })|_{y^{\ast }=0}$ . A good agreement with the experimental measurements is observed up to $x^{\ast }=105$  cm. The predicted maximum of $\unicode[STIX]{x1D70F}^{\ast }$ in the high-speed region and the minimum for the low-speed region occur at almost the same streamwise positions as in the measurements. It is rather remarkable that the nonlinear development of a primary instability (i.e. steady Görtler vortices) causes a skin friction increase that is a hallmark of laminar–turbulent transition. Lee & Liu (Reference Lee and Liu1992) found that nonlinear amplitude development, structures of Görtler vortices and the skin friction were quite well predicted by the solutions to the NBRE with an eigensolution as the initial condition. The comparisons in the present paper are more comprehensive and arguably more satisfactory in terms of the required input and accuracy. In our complete initial-boundary-value formulation, the only input to the calculation is the amplitude of FSVD. Use of eigensolutions as initial conditions could only be justified if there exists a discernible linear modal growth stage. That in turn requires the Görtler numbers to be reasonably large and FSVD intensity to be fairly low because otherwise, the linear modal stage is unattainable or bypassed as figure 4(a) indicates. The somewhat better accuracy (e.g. our figure 5 versus figure 4 in Lee & Liu (Reference Lee and Liu1992)) probably resulted from taking into account the continuing action of FSVD on the vortices in the boundary layer.

Figure 9. The nonlinear evolution of Görtler vortices as shown by contours of the streamwise velocity at different downstream locations for $G_{\unicode[STIX]{x1D6EC}}=1501$ and $R_{\unicode[STIX]{x1D6EC}}=1145$ . The increment of the streamwise velocity in the contour values is 0.1 $U_{\infty }$ .

Figure 9 presents the predicted evolution and structure of Görtler vortices by our computation using the experimental parameters of SB. Contour plots of the streamwise velocity, (4.3), at several streamwise locations are assembled to create a perspective view of the vortical structure. The continued action of the upwelling and downwelling motions leads to formation of mushroom-shaped patterns, which have been observed by many researchers (Hall Reference Hall1988; Lee & Liu Reference Lee and Liu1992; Li & Malik Reference Li and Malik1995; Ren & Fu Reference Ren and Fu2015). It is worth noting that streaks in the corresponding flat-plate boundary layer do not develop into a mushroom-shaped structure (Ricco et al. Reference Ricco, Luo and Wu2011).

4.2 Comparison with Ito (Reference Ito1985)

Figure 10. Comparison of the predicted streamwise velocity profiles with the experimental measurements (Ito Reference Ito1985). Solid and dashed lines: theoretical results at the peak and valley, respectively; squares and circles: experimental results at the peak and valley, respectively. The parameters are $G_{\unicode[STIX]{x1D6EC}}=460$ and $R_{\unicode[STIX]{x1D6EC}}=425$ .

Figure 11. The nonlinear evolution of Görtler vortices as shown by contours of the streamwise velocity in the range of $0.05<\hat{x}<1.0$ for $G_{\unicode[STIX]{x1D6EC}}=460$ and $R_{\unicode[STIX]{x1D6EC}}=425$ . The increment of the contour values is 0.1 $U_{\infty }$ .

The experiment of Ito (Reference Ito1985) was carried out for a boundary layer on a concave wall with a radius of curvature $r_{0}^{\ast }=1.0$  m, and the free-stream velocity $U_{\infty }=2.5~\text{ms}^{-1}$ . The spanwise wavelength of the Görtler vortices was $\unicode[STIX]{x1D706}_{z}\approx 1.6$  cm. On using these parameters, the corresponding Reynolds number $R_{\unicode[STIX]{x1D6EC}}=425$ and Görtler number $G_{\unicode[STIX]{x1D6EC}}=460$ . The parameters and the non-dimensional streamwise locations are shown in table 1. We choose $\unicode[STIX]{x1D716}=0.0022$ to fit the predicted amplitude with the experimental data at the first measurement position. The computational domain starts from $\hat{x}=0.001$ .

We calculated the streamwise velocity profiles at the peak and valley. The results are shown in figure 10, and compared with the experimental data of Ito (Reference Ito1985). The nonlinear effect is found to play a role from approximately $\hat{x}=0.35$ . A fairly good agreement is observed in both the linear and nonlinear development stages of the vortices. The predicted streamwise velocity profiles at the valley during nonlinear saturation at locations $\hat{x}=0.369$ , $0.462$ , $0.554$ and $0.646$ agree quite well with the measurements over almost the entire boundary layer. The agreement starts to deteriorate downstream of $\hat{x}=0.646$ , which suggests that secondary instability modes have become strong enough to influence the structure of saturated Görtler vortices.

Figure 11 presents the predicted evolution of Görtler vortices using Ito’s experimental parameters. As the values of $R_{\unicode[STIX]{x1D6EC}}$ and $G_{\unicode[STIX]{x1D6EC}}$ are smaller than those in SB: see table 1, a longer distance is required for the vortices to develop into the mushroom-shaped structure, the head of which is farther away from the wall, and the stem of the mushroom is longer.

4.3 Comparison with Mitsudharmadi et al. (Reference Mitsudharmadi, Winoto and Shah2004)

In the experiment of MWS, the surface has a radius of curvature $r_{0}^{\ast }=2$  m, and the free-stream velocity $U_{\infty }=3$ $\text{ms}^{-1}$ . Thirteen vertical wires were placed upstream of the leading edge of the test plate to generate Görtler vortices with a fixed spanwise wavelength, as opposed to naturally occurring ones in SB and Ito (Reference Ito1985). The spanwise distance between the wires is 15 mm, which sets the wavelength of the vortices. The corresponding Reynolds number $R_{\unicode[STIX]{x1D6EC}}=478$ and Görtler number $G_{\unicode[STIX]{x1D6EC}}=273$ . The presence of the vertical wires breaks isotropy, favouring components with longer length in the normal direction than that in the spanwise direction. For this reason, $k_{2}=0.6$ is adopted.

As a measure of the intensity of vortices, MWS defined the disturbance amplitude as

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D705}_{u}(\hat{x})=\underset{\unicode[STIX]{x1D702}}{\text{max}}(U_{D}(\unicode[STIX]{x1D702})-U_{U}(\unicode[STIX]{x1D702}))/(2U_{\infty }),\end{eqnarray}$$

where $U_{D}$ and $U_{U}$ denote the velocities at the downwash (peak) and upwash (valley) positions. The predicted amplitude will fit with the measured value at the first location if we take $\unicode[STIX]{x1D716}=0.0007$ , which corresponds to $Tu^{\dagger }=0.19\,\%$ , or $Tu=0.18\,\%$ if we take $\hat{u} _{1}^{\infty }=\hat{u} _{3}^{\infty }$ . The FST level in the test section is $0.35\,\%$ , five times as high as that in SB experiment. Again, this turbulence level apparently refers to unsteady fluctuations. Figure 12 compares $\unicode[STIX]{x1D705}_{u}$ with the measurement, and a good agreement is observed. It is evident that Görtler vortices are linear, i.e. undergo exponential amplification, in the region from $\hat{x}=0.1$ to 0.35, and then enter the nonlinear regime ranging from $\hat{x}=0.35$ to 0.75. Rapid growth of unsteady disturbances takes place from $\hat{x}\approx 0.72$ ( $x^{\ast }\approx 85$ cm), where Görtler vortices start meandering (MWS 2005). This is an indicator of secondary instability, which will be analysed in the next section.

Figure 13 displays a comparison of the streamwise velocity profiles with the experimental measurements. There is a fair qualitative agreement from $x^{\ast }=31$ cm to $90.4$ cm ( $\hat{x}=0.272$ to 0.792), but appreciable quantitative difference exists, even at the very first measurement point, where the vortices are weak and the base flow should be essentially Blasius. The difference suggests that a pressure gradient might be present. The discrepancy at $x^{\ast }=101.4$ cm ( $\hat{x}=0.888$ ) is most likely due to the onset and development of secondary instability.

4.4 Unsteady Görtler vortices

The response of the boundary layer over a curved plate to unsteady FSVD is now presented. We focus on the nonlinear evolution of unsteady Görtler vortices, and choose parameters $G_{\unicode[STIX]{x1D6EC}}=1501$ and $R_{\unicode[STIX]{x1D6EC}}=1145$ pertaining to the experiment of SB. The parameter $\hat{k}_{1}$ is related to the dimensional frequency $f^{\ast }$ via the relation

(4.7) $$\begin{eqnarray}\hat{k}_{1}=R_{\unicode[STIX]{x1D6EC}}k_{1}=2\unicode[STIX]{x03C0}R_{\unicode[STIX]{x1D6EC}}f^{\ast }\unicode[STIX]{x1D6EC}/U_{\infty }.\end{eqnarray}$$

Calculations were performed for three different frequencies $f^{\ast }=1$  Hz, 10 Hz and 20 Hz, corresponding to $\hat{k}_{1}=5.269$ , 52.691 and 104.382, respectively.

Figure 12. Comparison of the downstream development of the streamwise velocity fluctuation amplitude defined by (4.6): the line indicates the present computed results; the square symbols stand for the experimentally measured steady components (MWS 2005), whereas the circle symbols indicate the unsteady contributions.

Figure 13. Comparison of the streamwise velocity profiles with experimental measurements (Mitsudharmadi et al. Reference Mitsudharmadi, Winoto and Shah2005). Solid and dashed lines represent the theoretical predictions at the peak and valley, respectively; delta and square symbols represent the experimental results at the peak and valley, respectively. The parameters are $G_{\unicode[STIX]{x1D6EC}}=273$ , $R_{\unicode[STIX]{x1D6EC}}=478$ and $k_{2}=0.6$ .

Figure 14. Comparison of the downstream development of $\max u_{rms}$ for different $\hat{k}_{1}$ . The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501,R_{\unicode[STIX]{x1D6EC}}=1145$ and $\unicode[STIX]{x1D716}=0.0007$ .

Figure 15. Development of the fundamental mode $\max |\unicode[STIX]{x1D716}R_{\unicode[STIX]{x1D6EC}}2\hat{x}\hat{u} _{11}|$ and harmonic components $\max |\unicode[STIX]{x1D716}R_{\unicode[STIX]{x1D6EC}}2\hat{x}\hat{u} _{m,n}|$ : (a)  $\hat{k}_{1}=0$ ; (b)  $\hat{k}_{1}=5.269$ ; (c)  $\hat{k}_{1}=52.691$ ; (d)  $\hat{k}_{1}=104.382$ . The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501$ , $R_{\unicode[STIX]{x1D6EC}}=1145$ and $\unicode[STIX]{x1D716}=0.0007$ .

Figure 14 shows the development of $\max u_{rms}$ for these three frequencies with the steady case included for comparison. Vortices excited by lower-frequency disturbances acquire a larger amplitude prior to saturation. Interestingly, $\max u_{rms}$ saturates at almost the same level eventually. The saturation location shifts however downstream as the frequency increases. For the case of $\hat{k}_{1}=52.691$ , there appears a rather sudden change of the slope at $\hat{x}\approx 0.22$ . For the case of $\hat{k}_{1}=104.382$ , the first linear growth stage disappears and instead $\max u_{rms}$ decreases slightly after the transient growth. Starting from $\hat{x}=0.18$ , $\max u_{rms}$ grows exponentially at a fairly large rate until $\hat{x}\approx 0.32$ , where saturation starts.

In order to probe into the nature of nonlinear vortices, figure 15 shows the development of the maximum amplitudes of the fundamental and harmonic components for $\hat{k}_{1}=5.269,52.691$ and $104.382$ . For the lower-frequency case ( $\hat{k}_{1}=5.269$ ), the mean-flow distortion acquires a magnitude appreciably greater than that of the seeded fundamental modes $(1,\pm 1)$ ; a similar feature was observed in the case of the flat-plate boundary layer (Ricco et al. Reference Ricco, Luo and Wu2011). For the case of $\hat{k}_{1}=52.691$ , the amplitude of the mean-flow distortion $(0,0)$ is still larger than that of the fundamental modes, but the harmonic component (0,2) grows quickly from $\hat{x}\approx 0.23$ to overtake the fundamental modes $(1,\pm 1)$ . For the case $\hat{k}_{1}=104.382$ , the fundamental modes $(1,\pm 1)$ do not grow any more after the transient growth, instead the (0,2) mode amplifies exponentially. The mean-flow distortion (0,0) undergoes rapid growth again from about $\hat{x}=0.18$ to become the dominant component. It is worth noting that for both $\hat{k}_{1}=52.691$ and $104.382$ the rapid growth of the mean flow coincides with that of $\max u_{rms}$ shown in figure 14. In all the four cases, the (0,0) mode saturates at the highest level, and in the steady case the dominance of this mode is in accordance with the nonlinear theory of Hall & Lakin (Reference Hall and Lakin1988). For the steady and low-frequency cases, the fundamental modes acquire the second largest saturation amplitude, but once the frequency exceeds a critical value, it is the nonlinearly generated (0,2) mode that takes the second place.

Figure 16. Profiles of the instantaneous streamwise velocity $U$ at $z=0$ in different phases of time modulation: (a)  $\hat{k}_{1}=5.269$ ; (b)  $\hat{k}_{1}=52.691$ ; (c)  $\hat{k}_{1}=104.382$ . The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501,R_{\unicode[STIX]{x1D6EC}}=1145$ and $\unicode[STIX]{x1D716}=0.0007$ . The inserted window displays an enlarged view of the part in the dashed box.

Figure 17. Nonlinear evolution of Görtler vortices. (a,b)  $\hat{k}_{1}=5.269,\unicode[STIX]{x1D719}=25\unicode[STIX]{x03C0}/16$ (a) and $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ (b); (c,d)  $\hat{k}_{1}=52.691$ , $\unicode[STIX]{x1D719}=25\unicode[STIX]{x03C0}/16$ (c) and $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ (d); (e,f)  $\hat{k}_{1}=104.382,\unicode[STIX]{x1D719}=25\unicode[STIX]{x03C0}/16$ (e) and $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ (f). The increments of the contour levels is 0.1 $U_{0}$ . The parameters are $G_{\unicode[STIX]{x1D6EC}}=1501,R_{\unicode[STIX]{x1D6EC}}=1145$ and $\unicode[STIX]{x1D716}=0.0007$ .

Figure 16 shows the instantaneous streamwise velocity profiles at $z=0$ in different phases

(4.8) $$\begin{eqnarray}\unicode[STIX]{x1D719}\equiv \hat{k}_{1}\hat{\unicode[STIX]{x1D70F}},\end{eqnarray}$$

for three representative frequencies. For $\hat{k}_{1}=5.269$ , the velocity profile exhibits remarkable variation with the phase. The velocity may be larger or smaller than that of the Blasius flow. During certain time windows of the modulation, the instantaneous streamwise velocity profile becomes highly inflectional, suggesting that it may be inviscidly unstable. When $\hat{k}_{1}$ increases to $52.691$ , the temporal variation of the streamwise velocity becomes less strong, and the profile maintains a similar shape. The velocity near the wall hardly changes, and is persistently faster than the Blasius flow. However, a persistent large velocity deficit occurs at the outer part of the boundary layer; a similar feature was observed in the flat-plate case (Ricco et al. (Reference Ricco, Luo and Wu2011) and references therein). For the case $\hat{k}_{1}=104.382$ , the velocity profile changes moderately, i.e. Görtler vortices become almost stationary. This may be attributed to the dominance of the steady (0,0) and (0,2) modes. It is interesting to note that the profiles appear somewhat like the mean profile of fully developed turbulent boundary layer.

Figure 17 displays perspective views of the contours of the streamwise velocity at different downstream locations for $\hat{k}_{1}=5.269,52.691$ and $104.382$ . Two phases of time modulation $\unicode[STIX]{x1D719}=25\unicode[STIX]{x03C0}/16$ and $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ , are chosen. The perturbation eventually develops into a mushroom structure as in the steady case. For $\hat{k}_{1}=5.269$ , during a cycle of modulation, the mushrooms meander in the spanwise direction. For the case of $\hat{k}_{1}=52.691$ , the spanwise meandering is insignificant because the steady components, (0,0) and (0,2), dominate. Furthermore, since the (0,2) mode acquires a much larger amplitude than that of the fundamental, two mushrooms appear side by side within a fundamental spanwise wavelength. For the case of $\hat{k}_{1}=104.382$ , contour plots at the two phases are hardly distinguishable.

5.1 Formulation and methodology

Once Görtler vortices have acquired a sufficiently large amplitude, the distorted velocity profile $U(y,z;x,t)$ usually becomes inflectional in the streamwise and spanwise directions, and therefore is a harbinger of secondary instability (Saric Reference Saric1994), which amplifies high-frequency fluctuations and ultimately causes transition to turbulence. As was mentioned earlier, linear secondary instability analyses have been performed by several investigators (Hall & Horseman Reference Hall and Horseman1991; Yu & Liu Reference Yu and Liu1991; Li & Malik Reference Li and Malik1995). Previous calculations were carried out for Görtler vortices that developed spatially (or temporally) from somewhat arbitrary upstream (or initial) conditions, which suffices for the aim of establishing salient characteristics of the instability. In the present paper, secondary instability analyses will be conducted as part of the integrated approach. Our purpose is to add more quantitative information on the one hand, and build up the capacity for predicting transition in terms of given FST on the other.

Since the flow field of Görtler vortices varies with $\hat{x}$ and $\hat{\unicode[STIX]{x1D70F}}$ very slowly, the dependence on these two variables can be treated as being parametric when a short-wavelength (of order $\unicode[STIX]{x1D6FF}^{\ast }$ , with $\unicode[STIX]{x1D6FF}^{\ast }$ being the dimensional boundary-larger thickness) and high-frequency (of order $U_{\infty }/\unicode[STIX]{x1D6FF}^{\ast }$ ) instability is considered, just like the flat-plate case considered by Ricco et al. (Reference Ricco, Luo and Wu2011). With $\hat{x}$ and $\hat{\unicode[STIX]{x1D70F}}$ fixed, the base flow, comprising the Blasius flow and the Görtler vortices, depends only on $y$ and $z$ , and can be written as

(5.1) $$\begin{eqnarray}(\boldsymbol{U}_{G}(y,z),P_{G})=(U_{G}(y,z),V_{G}(y,z),W_{G}(y,z),P_{G}).\end{eqnarray}$$

The transverse velocities of the vortex base flow, $V_{G}$ and $W_{G}$ , will be neglected since they are very small, of $O(R_{\unicode[STIX]{x1D6EC}}^{-1})$ , and so $\boldsymbol{U}_{G}=(U,0,0)$ with $U$ being given by (4.3). When the base flow is perturbed by a disturbance of $O(\unicode[STIX]{x1D716}_{s})$ , the instantaneous flow is

(5.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{u}_{s}(x,y,z,t)=\boldsymbol{U}_{G}(y,z)+\unicode[STIX]{x1D716}_{s}\boldsymbol{u}_{s}^{\prime }(x,y,z,t),\\ p_{s}(x,y,z,t)=P_{G}(y,z)+\unicode[STIX]{x1D716}_{s}p_{s}^{\prime }(x,y,z,t),\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{u}_{s}^{\prime }=(u_{s}^{\prime },v_{s}^{\prime },w_{s}^{\prime })$ is the secondary perturbation velocity vector and $p_{s}^{\prime }$ the pressure. Substituting (5.2) into the N–S equations and neglecting the $O(\unicode[STIX]{x1D716}_{s}^{2})$ nonlinear terms, we obtain the linearized N–S equations,

(5.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}_{s}^{\prime }=-(\boldsymbol{U}_{G}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{s}^{\prime }-(\boldsymbol{u}_{s}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}_{G}-\unicode[STIX]{x1D735}p_{s}^{\prime }+\frac{1}{R_{\unicode[STIX]{x1D6EC}}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{s}^{\prime },\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}^{\prime }=0.\end{array}\right\}\end{eqnarray}$$

These equations, together with the homogeneous boundary conditions,

(5.4a,b ) $$\begin{eqnarray}\boldsymbol{u}_{s}^{\prime }=0\quad \text{at}~y=0,\quad (\boldsymbol{u}_{s}^{\prime },p_{s}^{\prime })\rightarrow 0\quad \text{as}~y\rightarrow \infty ,\end{eqnarray}$$

define a linear evolution operator $\ell$ , which evolves perturbations forward in time.

Given that the coefficients in (5.3) depend on $y$ and $z$ only, the solution for the secondary perturbation $\unicode[STIX]{x1D753}_{s}^{\prime }(x,y,z,t)\equiv (\boldsymbol{u}_{s}^{\prime },\,p_{s}^{\prime })$ can be sought of the modal form,

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D753}_{s}^{\prime }(x,y,z,t)=\unicode[STIX]{x1D753}_{s}(y,z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x-\unicode[STIX]{x1D714}t)}+\text{c.c}.,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D714}$ are the non-dimensional streamwise wavenumber and frequency of the secondary disturbance respectively. The shape function $\unicode[STIX]{x1D753}_{s}(y,z)$ is governed by a system of partial differential equations, ${\mathcal{L}}_{E}(\unicode[STIX]{x2202}_{y},\unicode[STIX]{x2202}_{z},U_{G};\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D714})\unicode[STIX]{x1D753}_{s}=0$ , which together with the homogeneous boundary conditions, $\boldsymbol{u}_{s}=0$ at $y=0$ and $(\boldsymbol{u}_{s},\,p_{s})\rightarrow 0$ as $y\rightarrow \infty$ , forms a so-called bi-global instability problem (Theofilis Reference Theofilis2011). In the temporal formulation, $\unicode[STIX]{x1D6FC}$ is real, whereas $\unicode[STIX]{x1D714}$ is complex valued and obtained as the eigenvalue. In the inviscid limit, the system for $\unicode[STIX]{x1D753}_{s}(y,z)$ can be reduced to a single equation for the pressure $p_{s}(y,z)$ (Hall & Horseman Reference Hall and Horseman1991),

(5.6) $$\begin{eqnarray}\left[\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}-\frac{2\unicode[STIX]{x2202}_{y}U}{U-\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}-\frac{2\unicode[STIX]{x2202}_{z}U}{U-\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}-\unicode[STIX]{x1D6FC}^{2}\right]p_{s}=0,\end{eqnarray}$$

where $U$ stands for the streamwise velocity (4.3) of the vortex flow. Despite that the secondary instability is essentially inviscid and that we have taken $R_{\unicode[STIX]{x1D6EC}}\gg 1$ in the analysis of the Görtler vortices, we will retain viscosity in the secondary instability calculations since this avoids the critical-layer singularity for neutral or nearly neutral modes, and furthermore may give more accurate results.

A discretization of the eigenvalue problem in two dimensions leads to a matrix, which is usually so large that it is difficult to compute its eigenvalues. To circumvent this obstacle, we use the Arnoldi method that makes use of time stepping. This method was developed by Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008), and was recently adapted to finite-difference setting by one of the present authors (Zhang & Luo Reference Zhang and Luo2015). The method is based on integrating the evolution system (5.3)–(5.4a,b ) from $t=0$ to $t=T$ with $T$ being a preset time. In implementation, we write $\unicode[STIX]{x1D719}_{s}^{\prime }(x,y,z,t)=\bar{\unicode[STIX]{x1D719}}(y,z,t)\text{e}^{\text{i}\unicode[STIX]{x1D6FC}x}+\text{c.c}.$ , where $\unicode[STIX]{x1D719}_{s}^{\prime }$ stands for any of $u_{s}^{\prime }$ , $v_{s}^{\prime }$ , $w_{s}^{\prime }$ and $p_{s}^{\prime }$ . Since the coefficients of the governing equations are periodic functions of $z$ , according to Floquet theory the solution for $\bar{\unicode[STIX]{x1D719}}$ takes the form,

(5.7) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}(y,z,t)=\text{e}^{\text{i}q\unicode[STIX]{x1D6FD}z}\mathop{\sum }_{n=-\infty }^{+\infty }\bar{\unicode[STIX]{x1D719}}^{(n)}(y,t)\text{e}^{\text{i}n\unicode[STIX]{x1D6FD}z},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}$ is the spanwise wavenumber and $0\leqslant q\leqslant 1/2$ . Only fundamental modes ( $q=0$ ) will be presented because the differences between the growth rates of the fundamental and subharmonic ( $q=1/2$ ) or detuned modes ( $0<q<1/2$ ) are less than $3.5\,\%$ (Ren & Fu Reference Ren and Fu2014). The discretization in the wall-normal direction is by a finite-difference scheme on variable grids, whose distribution is facilitated by a coordinate transformation from $y$ to $\tilde{\unicode[STIX]{x1D702}}$ (Zhang & Luo Reference Zhang and Luo2015), $y=y_{l}\tilde{\unicode[STIX]{x1D702}}/[1+b(1-\tilde{\unicode[STIX]{x1D702}}^{2})]$ with $b=(\sqrt{1+8k}-3)/4$ ; this maps the truncated physical domain $0\leqslant y\leqslant y_{l}$ to $0\leqslant \widetilde{\unicode[STIX]{x1D702}}\leqslant 1$ , where $k$ is the parameter controlling the mesh distribution. The interval (0,1) is divided into $N_{y}$ equally spaced subintervals, and so $\unicode[STIX]{x0394}\widetilde{\unicode[STIX]{x1D702}}=1/N_{y}$ . Typically, we use 300 points in the normal direction, and 64 terms in the Fourier series, but for modes concentrating in the wall region, refined grids are needed there in order to resolve the eigenfunctions.

Let $\bar{\unicode[STIX]{x1D753}}$ denote the discretized solution vector. The integration from $t=0$ to $t=T$ , starting with an initial perturbation $\bar{\unicode[STIX]{x1D753}}_{0}$ , may be viewed as a mapping,

(5.8) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D753}}(T)=\ell (T)\bar{\unicode[STIX]{x1D753}}_{0},\end{eqnarray}$$

by the evolution operator $\ell$ . Repeating this, one obtains $\ell ^{2}\bar{\unicode[STIX]{x1D753}}_{0},\ldots ,\ell ^{k}\bar{\unicode[STIX]{x1D753}}_{0}$ , which can be arranged to form the normalized Krylov sequence,

(5.9) $$\begin{eqnarray}{\mathcal{T}}_{k+1}=[\bar{\unicode[STIX]{x1D753}}_{0}/a_{0},\ell \bar{\unicode[STIX]{x1D753}}_{0}/a_{1},\ldots ,\ell ^{k}\bar{\unicode[STIX]{x1D753}}_{0}/a_{k}],\end{eqnarray}$$

where $a_{k}$ are the normalization factors to make each column a unit vector. A simple fundamental relation holds between ${\mathcal{T}}_{k+1}$ and $\ell (T){\mathcal{T}}_{k}$ . Using this relation and carrying out QR decompositions of ${\mathcal{T}}_{k+1}$ and ${\mathcal{T}}_{k}$ , a $k\times k$ Hessenberg matrix $\unicode[STIX]{x1D643}_{k}$ can be constructed. Since $k$ is rather small, eigenvalues of $\unicode[STIX]{x1D643}_{k}$ , $\unicode[STIX]{x1D706}_{j}$ ( $j=1,2,\ldots ,k$ ), and the corresponding eigenvectors, can be computed at little cost. The eigenvalues for the original problem, $\unicode[STIX]{x1D714}$ , are related to $\unicode[STIX]{x1D706}_{j}$ by $\unicode[STIX]{x1D714}=(\text{i}/T)\ln \unicode[STIX]{x1D706}_{j}$ , and the eigenfunctions can be constructed by following the algorithm described in Barkley et al. (Reference Barkley, Blackburn and Sherwin2008).

5.2 Secondary instability of steady Görtler vortices

5.2.1 Secondary instability calculations pertaining to the SB (1987) experiment

Although multiple unstable eigenmodes co-exist for each Görtler vortex flow profile, we focus on the most unstable ones. Both the computational results of the last section and the measurements of SB show that at $\hat{x}=0.227$ ( $x^{\ast }=95$ cm) Görtler vortices acquire their maximum amplitude. The experiment indicates that unsteady disturbances amplify abruptly from this position. The abrupt amplification is most likely due to the development of secondary instability modes, triggered by high-frequency components in FSVD, but the precise receptivity mechanism remains a complete mystery at the present. We select nine streamwise locations (listed in table 1), at which secondary stability analysis is to be performed.

Figure 18. Characteristics of secondary instability of Görtler vortices in the experiment of SB: (a) the temporal growth rate and (b) the phase speed versus the streamwise wavenumber $\unicode[STIX]{x1D6FC}$ at $x^{\ast }=$ 80, 85, 87 and 90 cm, corresponding to $\hat{x}=0.191,0.203$ , $0.209$ and $0.215$ , respectively.

Figure 19. Eigenfunctions (absolute value, black line) of secondary unstable modes, shown by contours of the streamwise velocity at $\hat{x}=0.215$ . (a) The most unstable mode with $\unicode[STIX]{x1D6FC}=0.8$ (odd mode I); (b) the most unstable mode with $\unicode[STIX]{x1D6FC}=1.5$ (even mode I). The grey contours represent the streamwise velocity of the vortex base flow. Nine levels are specified, ranging from 0.1 to 0.9.

Figure 18 plots the growth rates and phase speeds at $\hat{x}=0.191,0.203,0.209$ and $0.215$ . Figure 18(a) shows that at $\hat{x}=0.191$ , only unstable modes are found for $\unicode[STIX]{x1D6FC}<1$ , and these are sinuous modes, or odd modes I, as referred to by Li & Malik (Reference Li and Malik1995). Another family of modes with larger $\unicode[STIX]{x1D6FC}$ , referred to as even modes I by Li & Malik (Reference Li and Malik1995), start to emerge from $\hat{x}=0.203$ , and their growth rates increase. However, odd modes I remain dominant up to $\hat{x}=0.215$ . The maximum growth rates at the first three streamwise locations are attained for streamwise wavenumbers approximately between 0.6 and 0.8, corresponding to wavelengths between 3.8 and 2.9 cm. The growth rates of even modes I become almost equal to those of the odd modes I when $\hat{x}=0.215$ . The result in the figure also indicates that the band of unstable modes becomes broader at downstream positions. Figure 18(b) shows phase speeds of secondary instability modes. The phase speeds of sinuous modes decrease rather rapidly with $\unicode[STIX]{x1D6FC}$ , whereas those of even modes I exhibit a weaker dependence on $\unicode[STIX]{x1D6FC}$ . The variation of the growth rates with $\unicode[STIX]{x1D6FC}$ and streamwise location, is similar to the inviscid predictions of Li & Malik (Reference Li and Malik1995).

Figure 19 displays the eigenfunctions of typical sinuous and varicose modes in terms of their streamwise velocity. The eigenfunction is normalized such that its maximum value is unity. As figure 19(a) indicates, an odd mode I has two pairs of dominant peaks, one pair is located at the outer part of the boundary layer, while the other is well within the boundary layer. The modal shape is similar to the inviscid predictions of Hall & Horseman (Reference Hall and Horseman1991) and Li & Malik (Reference Li and Malik1995). A notable difference is that the inner peaks in Hall & Horseman (Reference Hall and Horseman1991) are farther out from the valley centre, whereas those in our calculation and in Li & Malik (Reference Li and Malik1995) are nearer to the valley centre, resembling more closely the measurement shown in figure 16 of SB. All these peaks appear to reside in the region where $U$ varies quickly with $z$ . Figure 19(b) shows that an even mode I has three dominant peaks, which are in the region where $U$ varies quickly with $y$ .

Figure 20. Characteristics of secondary instability: (a) the temporal growth rate and (b) the phase speed versus the streamwise spanwise wavenumber $\unicode[STIX]{x1D6FC}$ at $x^{\ast }=90$ , 95, 100 and 105 cm corresponding to $\hat{x}=0.215,0.227,0.238$ and $0.250$ , respectively.

Figure 21. The characteristics of a sinuous mode II with a wavenumber $\unicode[STIX]{x1D6FC}=2.9$ at $\hat{x}=0.227$ ( $x^{\ast }=95$  cm). (a) Contours of the streamwise velocity eigenfunction (absolute value, solid line) and the base flow (dashed lines). Nine contour levels are specified ranging from 0.1 to 0.9. (b) The variation of growth rate (solid line) and frequency (dashed line) with $R_{\unicode[STIX]{x1D6EC}}$ ; the circle indicates the position $R_{\unicode[STIX]{x1D6EC}}=1145$ .

Figure 20 shows the growth rates and phase speeds at $\hat{x}=0.215$ , $0.227$ , $0.238$ and $0.250$ . Both even modes I and odd modes I are found at these four streamwise positions further downstream. The stability characteristics changes at $\hat{x}=0.215$ in that downstream of this position, there emerge odd modes II, which have larger wavenumbers. The existence of this family of modes was noted only recently by Ren & Fu (Reference Ren and Fu2015). From $\hat{x}=0.227$ , for $\unicode[STIX]{x1D6FC}$ less than approximately 1.0, odd modes I are dominant at all these streamwise positions, but for $\unicode[STIX]{x1D6FC}>1.0$ , even modes I dominate until $\unicode[STIX]{x1D6FC}$ reaches a certain value, where odd modes II take over. The growth rates of secondary instability modes increase until $\hat{x}=0.250$ , after which they drop as Görtler vortices attenuate, but interestingly the spectral range of unstable modes continues to expand. Figure 20(b) shows that the phase speeds of odd modes II are lower than those of other modes at every streamwise position. Hence odd modes II tend to concentrate in the wall region and may be difficult to detect, which might be the reason why they were not discovered by SB (1987).

The eigenfunction of a typical odd mode II is shown in figure 21(a). The streamwise velocity $|u_{s}|$ of the odd mode II at $x^{\ast }=95$ cm with a wavenumber $\unicode[STIX]{x1D6FC}=2.9$ is plotted. A question arises as to whether odd modes II are of viscous nature or essentially inviscid. In order to answer this question, we calculated the growth rates $\unicode[STIX]{x1D714}_{i}$ and frequencies $\unicode[STIX]{x1D714}_{r}$ of odd modes II for very high Reynolds numbers, which represent the inviscid asymptotic limit. As is shown in figure 21(b), the growth rate increases monotonically with $R_{\unicode[STIX]{x1D6EC}}$ and finally approaches a constant, and the frequency $\unicode[STIX]{x1D714}_{r}$ becomes independent of $R_{\unicode[STIX]{x1D6EC}}$ too, implying that odd modes II are inviscid in character.

Figure 22. (a,b) The spatial growth rates of secondary instability versus the frequency at different streamwise locations.

When a temporal secondary mode ( $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D714}$ ) is found, the spatial growth rate can be obtained from the temporal rate by using Gaster transformation (Gaster 1962), according to which the complex wavenumber $\tilde{\unicode[STIX]{x1D6FC}}$ is related to $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D714}$ via

(5.10) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6FC}}=\tilde{\unicode[STIX]{x1D6FC}_{r}}+\text{i}\tilde{\unicode[STIX]{x1D6FC}_{i}}=\left(\unicode[STIX]{x1D6FC}-\frac{\unicode[STIX]{x1D714}_{i}}{|C_{g}|}\sin \unicode[STIX]{x1D712}\right)-\text{i}\left(\frac{\unicode[STIX]{x1D714}_{i}}{|C_{g}|}\cos \unicode[STIX]{x1D712}\right),\end{eqnarray}$$

where $C_{g}=\unicode[STIX]{x2202}\unicode[STIX]{x1D714}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}\equiv |C_{g}|\text{e}^{\text{i}\unicode[STIX]{x1D712}}$ is the (complex-valued) group velocity with $\unicode[STIX]{x1D712}$ being the phase of $C_{g}$ . Figure 22 displays the spatial growth rates so obtained at several locations in the range of $x^{\ast }=80\sim 105$ cm. Lower frequency unstable waves ( $\unicode[STIX]{x1D714}_{r}<0.7$ ) grow first but at smaller rates. The dimensional frequency $f^{\ast }$ reported in experiments is related to $\unicode[STIX]{x1D714}_{r}$ by

(5.11) $$\begin{eqnarray}f^{\ast }=\unicode[STIX]{x1D714}_{r}U_{\infty }/(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6EC}).\end{eqnarray}$$

Upstream of $\hat{x}=0.215$ , the most dangerous sinuous modes are centred at $\unicode[STIX]{x1D714}_{r}=0.6$ , corresponding to $f^{\ast }=128$  Hz. In the experiment of SB, the most dominant wave has a frequency of 130 Hz, which is very close to our theoretical prediction. In addition, the iso-contours of the streamwise velocity fluctuations at $\hat{x}=0.191,0.215$ and $0.238$ in figure 16 of SB, exhibit a structure very similar to that of odd modes I. Comparisons were made previously by Sabry & Liu (Reference Sabry and Liu1991) and Yu & Liu (Reference Yu and Liu1991), but their instability analysis was performed for temporally evolving vortices. For spatially evolving vortices, the inviscid calculations (Hall & Horseman Reference Hall and Horseman1991; Li & Malik Reference Li and Malik1995) gave the instability characteristics (including the range of frequency) broadly in agreement with experiments. The amplification of secondary instability modes may explain the abrupt change of the intensity of unsteady disturbances, shown in figure 6.

Figure 22 also shows that the growth rate of the most unstable higher-frequency even mode ( $\unicode[STIX]{x1D714}=0.8$ ) is comparable with that of the lower-frequency one ( $\unicode[STIX]{x1D714}=0.4$ ) at $\hat{x}=0.215$ . The frequency of the former is about twice that of the latter, and will be even higher at locations farther downstream. The band of unstable modes becomes broader downstream despite the fact that the intensity of Görtler vortices decreases, implying that not only the amplitude but also the spatial structure of vortices is relevant to the instability characteristics.

Figure 23. The local spatial growth rates of modes each with a fixed frequency, and comparison with the experimental data of SB (1987). The results are shown for four different frequencies, $\unicode[STIX]{x1D714}=0.3$ (65 Hz), $0.5$ (107 Hz), $0.598$ (130 Hz) and $0.8$ (174 Hz).

We also traced the downstream development of several modes, each with a fixed frequency. The variation of their growth rates with the streamwise distance is displayed in figure 23. Obviously, the mode with $\unicode[STIX]{x1D714}=0.598$ ( $f^{\ast }=130$  Hz) starts to amplify from $\hat{x}=0.190$ ( $x^{\ast }=80$  cm). Its maximum growth occurs in the range of $\hat{x}=0.225\sim 0.250$ ( $x^{\ast }=94\sim 105$  cm), which is same as the experimental observation. Also displayed in figure 23 is the growth rate of the r.m.s. of unsteady perturbations, calculated by

(5.12) $$\begin{eqnarray}\unicode[STIX]{x1D70E}=\frac{\unicode[STIX]{x1D6EC}}{u_{rms}^{\prime }}\frac{\unicode[STIX]{x2202}u_{rms}^{\prime }}{\unicode[STIX]{x2202}x^{\ast }},\end{eqnarray}$$

where $u_{rms}^{\prime }$ denotes the experimentally measured r.m.s. of total unsteady perturbations. Since no experimental data on growth rates of secondary instability are available, we take $\unicode[STIX]{x1D70E}$ as their proxy. The predicted spatial growth rates $-\unicode[STIX]{x1D6FC}_{i}$ are more or less of the same order of magnitude as the measured $\unicode[STIX]{x1D70E}$ , but the former are several times larger than the latter. This quantitative discrepancy is not surprising because $-\unicode[STIX]{x1D6FC}_{i}$ and $\unicode[STIX]{x1D70E}$ are, though related, not the same quantity. In experiment, rather than following the development of a well-defined mode, the growth rate $\unicode[STIX]{x1D70E}$ of unsteady disturbances was computed from $u_{rms}$ , which includes a broad band of components. The unstable mode is smeared, and thus $\unicode[STIX]{x1D70E}$ , which measures the amplification rate of overall unsteady disturbances, is expected to be smaller than $-\unicode[STIX]{x1D6FC}_{i}$ .

The results in figures 22 and 23 indicate that higher-frequency even modes start to amplify later but at larger rates. Their presence and intensity in the boundary layer depend, like all modes, on the content of high-frequency FSVD and receptivity. If triggered, on the theoretical ground they are expected eventually to dominate farther downstream. We note that such even modes, which manifest themselves as horseshoe structures, were indeed observed by Aihara & Koyama (Reference Aihara and Koyama1981) and Ito (Reference Ito1985). However, in the experiment of SB even modes appeared rather infrequently, probably because they have not yet gained significant amplitudes in the region of measurement.

5.2.2 Secondary instability calculations pertaining to MWS (2005) experiment

MWS (2005) studied the onset of secondary instability of forced wavelength Görtler vortices by experiment. Through a spectral analysis of velocity fluctuations, they found that the peak frequency is approximately 150 Hz, and the spanwise wavelength is comparable with the spacing of the primary Görtler vortices, implying that the fundamental secondary instability is the dominant mode in their experiment. Secondary instability calculations are performed for Görtler vortices pertaining to this condition.

Figure 24. The characteristics of secondary instability of Görtler vortices in the experiment of MWS (2005): (a) the temporal growth rates, and (b) the phase speeds versus the streamwise wavenumber $\unicode[STIX]{x1D6FC}$ at different streamwise locations.

Figure 25. The spatial growth rate versus frequency at different streamwise locations.

Figure 26. The local spatial growth rates of modes each with a fixed frequency and comparison with the experimental data of MWS (2005). Results are shown for three unstable waves with $\unicode[STIX]{x1D714}=0.5$ (100 Hz), $0.75$ (150 Hz) and $1.5$ (300 Hz).

Figure 24(a) displays the temporal growth rate for a range of the streamwise wavenumbers at different downstream locations. In the upstream region up to $x^{\ast }\approx 60$  cm, only odd modes I exist. Three families of modes, odd modes I, even modes I and odd modes II, are observed at $x^{\ast }=65$  cm. The growth rates of even modes increase until $x^{\ast }=70$  cm, after which they decrease. From $x^{\ast }=70$  cm, the growth rates of odd modes II increase to attain their respective maximum values, while those of odd modes I decrease slightly. The range of most unstable even modes I shrinks. At $x^{\ast }=95.2$  cm, odd modes II become most unstable. Figure 24(b) shows the phase speeds of the most unstable modes identified in different wavenumber ranges. The temporal growth rates are again converted into spatial ones using Gaster transformation. Figure 25 plots the spatial growth rates so obtained versus the frequency at different streamwise locations. The dominant frequencies of the odd mode I, even mode I and odd mode II are found to be 0.5 (100 Hz), 0.75 (150 Hz) and 1.5 (260 Hz), respectively. Obviously, large temporal growth rates translate broadly into large spatial amplification. Since strong spatial growth rates persisting over a streamwise region imply substantial accumulated amplification (see below), which determines which mode at a given location would dominate, it may be inferred that as vortices develop downstream the dominant disturbances are odd modes I first, then even modes follow, and odd modes II will dominate farther downstream. That odd modes I grow first is consistent with the experiment observation. MWS found that the fundamental frequency was in the range of 120–180 Hz, which is in reasonable agreement with our theoretical result. Given the eigenfunction of an odd mode II, fluctuations are expected to concentrate in the region near the ‘stem’ of the mushroom prior to turbulence. The turbulent intensity contours shown in figure 10 of Mitsudharmadi et al. (Reference Mitsudharmadi, Winoto and Shah2004) are consistent with this expectation, and may be taken as evidence for the existence of odd modes II.

Figure 26 shows the variation with the streamwise distance of the spatial growth rates of secondary instability modes with fixed frequencies $\unicode[STIX]{x1D714}=0.5$ (100 Hz), $0.75$ (150 Hz) and $1.5$ (300 Hz). As for the case of SB, lower-frequency waves grow first and expected to be observed upstream, but higher-frequency waves are likely to become significant too due to their accumulated amplification, over the regions $x^{\ast }>60$ and $x^{\ast }>65$ respectively for the two modes shown. Again, the predicted growth rates of secondary instability are compared with the measured $\unicode[STIX]{x1D70E}$ , and they are found to be of the same order of magnitude. However, significant difference exists, a possible reason for which was given earlier.

5.3 Secondary instability of unsteady Görtler vortices

During the evolution of unsteady Görtler vortices, the distorted instantaneous velocity profiles appear, in certain phases, inflectional in the streamwise and spanwise directions, and thus may become inviscidly unstable. This expectation is confirmed by our secondary instability analysis.

Figure 27. Instability characteristics of unsteady Görtler vortices with $\hat{k}_{1}=5.291$ at $\hat{x}=0.215$ ( $x^{\ast }=90$  cm) (a,b) and $\hat{x}=0.238$ ( $x^{\ast }=100$  cm) (c,d) the temporal growth rate $\unicode[STIX]{x1D714}_{i}$ and the phase speed $C_{r}$ versus the wavenumber $\unicode[STIX]{x1D6FC}$ at seven instants.

Figure 27(a,b) shows the instability characteristics at $\hat{x}=0.215$ ( $x^{\ast }=90$  cm), where unsteady Görtler vortices with frequency $\hat{k}_{1}=5.29$ (2 Hz) reach their maximum amplitude. Calculations were performed at seven instants corresponding to $\unicode[STIX]{x1D719}=13\unicode[STIX]{x03C0}/16$ , $\unicode[STIX]{x03C0},20\unicode[STIX]{x03C0}/16,22\unicode[STIX]{x03C0}/16$ , $25\unicode[STIX]{x03C0}/16,30\unicode[STIX]{x03C0}/16$ and $31\unicode[STIX]{x03C0}/16$ , the streamwise velocity profiles at which are displayed in figure 16. Within a short time interval around $\unicode[STIX]{x1D719}=25\unicode[STIX]{x03C0}/16$ , only a small band of weakly growing modes exists with $\unicode[STIX]{x1D6FC}$ near 0.9. Vigorous growth occurs around $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ and $31\unicode[STIX]{x03C0}/16$ . The secondary instability is thus intermittent in time, similar to the flat-plate case studied by Wu & Choudhari (Reference Wu and Choudhari2003) and Ricco et al. (Reference Ricco, Luo and Wu2011). The maximum growth rate is more than twice that of the steady case shown in figure 20 at the same streamwise location. Figure 27(b) plots the phase speeds, which fall in the range of $0.5{-}0.7U_{\infty }$ . Figure 27(c,d) presents the growth rates and phase speeds at $\hat{x}=0.238$ ( $x^{\ast }=100$  cm). Both the growth rate and the phase speed are slightly larger in all phases than those at $\hat{x}=0.215$ . The temporal growth rates are merely suggestive of what kind of disturbances are likely to emerge. As in the steady case, the temporal growth rates can be converted to spatial ones, and the actual presence and intensity of the secondary disturbance depend on receptivity. However, for unsteady Görtler vortices the secondary disturbance is modulated both in time and space, tracing its development is more complicated than integrating the local spatial growth rate over the streamwise distance, and is not pursued in the present paper.

Figure 28. Contours of the eigenfunction (the streamwise velocity) at $\hat{x}=0.238$ ( $x^{\ast }=100$  cm) for $\hat{k}_{1}=5.691$ : (a–c) $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ ; (a)  $\unicode[STIX]{x1D6FC}=1.2$ , (b)  $\unicode[STIX]{x1D6FC}=2.3$ , (c)  $\unicode[STIX]{x1D6FC}=3.7$ ; (d–f) $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ; (d)  $\unicode[STIX]{x1D6FC}=0.9$ , (e)  $\unicode[STIX]{x1D6FC}=2.1$ , (f)  $\unicode[STIX]{x1D6FC}=3.8$ .

Figure 28 plots the eigenfunctions of secondary instability at $\hat{x}=0.238$ , and contours of the instantaneous streamwise velocity of Görtler vortices are also included. Similar to the steady case, odd modes I, even modes I and odd modes II exist for small, intermediate and large $\unicode[STIX]{x1D6FC}$ . Note that for this frequency, unsteady Görtler vortices meander in the spanwise direction, and so do the eigenfunctions of secondary instability.

Next we consider Görtler vortices with a moderate frequency $\hat{k}_{1}=52.916$ ( $f^{\ast }=10$  Hz). Since saturation occurs later than in the steady and lower-frequency cases, secondary instability calculations are carried at positions somewhat farther downstream, at $\hat{x}=0.288$ and $0.336$ ( $x^{\ast }=120$  cm and $140$  cm). The results are shown in figure 29. The variation of growth rates and phase speeds of dominant odd modes I becomes moderate, i.e. the intermittency of secondary instability weakens. The growth rates of odd modes II, centred at the second peak of these curves, still vary considerably with the phase at $\hat{x}=0.288$ , but their variation becomes weaker too at $\hat{x}=0.336$ even though their phase speeds exhibit strong dependence on the phase (figure 29 d).

Figure 29. Instability characteristics of unsteady Görtler vortices with $\hat{k}_{1}=52.916$ at $\hat{x}=0.288$ ( $x^{\ast }=120$  cm) (a,b) and $\hat{x}=0.336$ ( $x^{\ast }=140$  cm) (c,d) the temporal growth rate and the phase speed $C_{r}$ versus $\unicode[STIX]{x1D6FC}$ .

Figure 30 displays the eigenfunctions of secondary instability at $\hat{x}=0.336$ ( $x^{\ast }=140$ cm) for $\hat{k}_{1}=52.691$ . At each of the two instants, eigenfunctions of an odd mode I, an even mode I and an odd mode II are shown. The spanwise wavelength of the odd mode I is half of the fundamental period of the vortices, indicating that odd modes are controlled primarily by the strong (0,2) component. They concentrate in the shoulder and stem regions of Görtler vortices. The even mode I and odd mode II have the same fundamental periodicity. They are intermittent, and reside on the top and in the stem region of the mushroom respectively.

Figure 30. Contours of the eigenfunctions (the streamwise velocity) at $\hat{x}=0.336$ ( $x^{\ast }=140$ cm) for $\hat{k}_{1}=52.691$ : (a–c) $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ ; (a)  $\unicode[STIX]{x1D6FC}=1.3$ , (b)  $\unicode[STIX]{x1D6FC}=2.0$ , (c)  $\unicode[STIX]{x1D6FC}=3.1$ ; (d–f) $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ; (d)  $\unicode[STIX]{x1D6FC}=1.2$ , (e)  $\unicode[STIX]{x1D6FC}=2.1$ , (f)  $\unicode[STIX]{x1D6FC}=3.0$ .

Figure 31. Instability characteristics of unsteady Görtler vortices with $\hat{k}_{1}=104.382$ at $\hat{x}=0.360$ ( $x^{\ast }=150$  cm) (a,b) and $\hat{x}=0.381$ ( $x^{\ast }=160$  cm) (c,d) the temporal growth rate and the phase speed $C_{r}$ versus $\unicode[STIX]{x1D6FC}$ at different phases.

Figure 31 shows the secondary instability characteristics of Görtler vortices with $\hat{k}_{1}=104.382$ ( $f^{\ast }=20$  Hz) at $\hat{x}=0.360$ and $0.381$ ( $x^{\ast }=150$  cm and $160$  cm). For such a high frequency, the intermittence of secondary instability almost disappears. Both the growth rates and phase speeds vary very little with the phase at each streamwise location. The maximum growth rates of odd and even modes I at $\hat{x}=0.381$ become larger than those at $\hat{x}=0.360$ , but the growth rates of odd modes II decrease. Contours of the eigenfunction (the streamwise velocity) at $\hat{x}=0.381$ ( $x^{\ast }=160$  cm) are shown in figure 32. For each of three types of modes, the eigenfunctions at the two phases appear almost identical (steady). Their spanwise wavelength is half of the fundamental period. We may conclude that the secondary instability is completely controlled by the (0,2) component.

Figure 32. Contours of the eigenfunction (the streamwise velocity) at $\hat{x}=0.381$ ( $x^{\ast }=160$  cm) for $\hat{k}_{1}=104.382$ : (a–c) $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ; (a) $\unicode[STIX]{x1D6FC}=1.3$ , (b)  $\unicode[STIX]{x1D6FC}=1.9$ , (c)  $\unicode[STIX]{x1D6FC}=3.0$ ; (d–f) $\unicode[STIX]{x1D719}=31\unicode[STIX]{x03C0}/16$ ; (d)  $\unicode[STIX]{x1D6FC}=1.3$ , (e)  $\unicode[STIX]{x1D6FC}=2.0$ , (f)  $\unicode[STIX]{x1D6FC}=3.0$ .