2.1. Curved channel problem

Consider an incompressible fluid with kinematic viscosity $\nu ^{*}$ flowing through a non-uniformly curved channel of dimensional width $2d^{*}$ . The channel curvature is specified by $a^{*}$ , the dimensional local radius of curvature of the channel centreline. It is assumed that the channel is not curved in the spanwise direction. We employ an orthogonal curvilinear coordinate system $(x^{*}, y^{*}, z^{*})$ aligned with the geometry of the channel. Here, $y^{*}$ and $z^{*}$ measure distances along two orthogonal directions within the cross-section, while $x^{*}$ is defined as the arc-length coordinate along the channel centreline, starting from the inlet. In other words, we adopt body-fitted coordinates (Slattery Reference Slattery1999) attached to the channel centreline. As shown in figure 1, the $y^{*}$ -axis is taken normal to the wall, while the $z^{*}$ -axis is aligned in the spanwise direction. To avoid confusion, Cartesian coordinates that share the origin with the body-fitted coordinates are denoted by $(\check {x}^{*},\check {y}^{*},\check {z}^{*})$ .

Figure 1. Schematic of the curved channel flow problem. Length is non-dimensionalised using the half-width of the channel. The body-fitted coordinates $(x,y,z)$ are attached to the channel centreline. For a given radius of curvature $a$ , the centreline has a unique parametrisation in Cartesian coordinates $(\check {x},\check {y},\check {z})$ . In §§ 4.1 and 4.2, the inlet conditions are travelling-wave ECS of plane Poiseuille flow (see § 3.2). In § 4.3, the inlet and outlet are reversed, and the inlet condition is the Dean vortex solution obtained in § 4.2.

Using the half-width $d^{*}$ as the length scale, the bulk velocity $U_b^{*}$ as the velocity scale and the fluid density as the density scale, the governing equations in the non-dimensional coordinates $(x,y,z)=(x^{*},y^{*},z^{*})/d^{*}$ are written as

(2.1) \begin{eqnarray} (\partial _t +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}=-\boldsymbol{\nabla \!}p+R^{-1}\boldsymbol{\nabla} ^2 \boldsymbol{u},\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{eqnarray}

where $R=U_b^{*}d^{*}/\nu ^{*}$ is the bulk Reynolds number. The pressure gradient $\boldsymbol{\nabla \!}p$ includes a component that ensures the non-dimensional bulk velocity remains unity. The velocity field can be expressed as $\boldsymbol{u}=u\boldsymbol{e}_x+v\boldsymbol{e}_y+w\boldsymbol{e}_z$ , where $\boldsymbol{e}_x$ , $\boldsymbol{e}_y$ and $\boldsymbol{e}_z$ are the unit vectors in the $x$ -, $y$ - and $z$ -directions, respectively. The component form of the governing equations in body-fitted coordinates is provided in Appendix A.

Throughout this paper, we consider the weak curvature regime introduced by Dean (Reference Dean1928), in which $a^{*}$ is assumed to be $O(R^{2}d^{*})$ . Following convention, we introduce the Dean number, $De$ , such that the non-dimensional radius of curvature is given by

(2.2) \begin{eqnarray} a=a^{*}/d^{*}=8R^2/De^2. \end{eqnarray}

Similar to Hall’s analysis of Görtler vortices (Hall Reference Hall1983, Reference Hall1988), we assume that $a$ varies slowly with $R^{-1}x$ .

The channel walls are located at $y=\pm 1$ , where the no-slip boundary conditions are imposed. In the $z$ direction, periodicity with period $2\pi /\beta$ is assumed, where $\beta$ denotes the spanwise wavenumber. To fully define the system, appropriate boundary conditions must also be specified at the upstream and downstream ends of the computational domain. Our interest lies in time-periodic solutions (i.e. global ECS) of this system.

The need to specify upstream and downstream boundary conditions poses a difficulty not encountered in ECS computations with periodic domains. These conditions must be time periodic; otherwise, a time-periodic solution of the system cannot exist. In the scenario we consider, the curvature is assumed to be constant sufficiently far upstream of the inlet. This assumption allows us to use results obtained under periodic boundary conditions there. For example, when the curvature is zero upstream, the flow there reduces to plane Poiseuille flow, where travelling-wave-type ECS are well defined. In this setting, the inlet boundary condition naturally becomes time periodic. As noted in § 1, generating travelling waves in experiments is not straightforward. Nevertheless, given the established importance of ECS to the dynamics of plane Poiseuille flow, this choice of inlet condition is reasonable. On the other hand, no boundary conditions are required at the outlet, since we will solve a parabolic reduced form of the governing equations that can be integrated in the streamwise direction.

Section 2.2 briefly reviews the derivation of the BRE, which represents the simplest form of parabolised equations. We shall see that the BRE are limited to capturing fluid motions that evolve on the viscous time scale, $t^{*}=O(d^{*2}/\nu ^{*})$ , and over the corresponding diffusive streamwise length scale, $x^{*}=O(U_b^{*}d^{*2}/\nu ^{*})$ .

In contrast, self-sustained travelling waves in plane Poiseuille flow (i.e. local ECS at the inlet) typically exhibit dynamics on the convective spatio-temporal scale, $t^{*}=O(d^{*}/U_b^{*})=O(R^{-1}d^{*2}/\nu ^{*})$ and $x^{*}=O(d^{*})=O(R^{-1}U_b^{*}d^{*2}/\nu ^{*})$ . Therefore, in § 2.3, we propose a new computational approach that reconciles the slow scales captured by the BRE with the fast scales associated with the SSP-VWI. This novel approach combines the idea of parabolisation with the computation of ECS, and is hereafter referred to as the parabolised coherent structure (PCS) approach. Our computational method, presented in § 3, and all results in § 4 are based on the PCS system, which is derived under the four assumptions described in § 2.3.

While the BRE are directly derived via large-Reynolds-number asymptotic analysis, the PCS system is not. Nevertheless, as we shall see in § 2.4 and Appendix A, the PCS system retains all terms that are important in the asymptotic analysis combining fast- and slow-scale dynamics. In other words, the PCS equations include terms smaller than leading order. These terms do not interfere with capturing the leading-order fluid motion, yet they make numerical computations much simpler than solving the leading-order system directly. Distinguishing leading- and higher-order terms requires rather intricate analysis, because the dominant balance differs inside and outside the thin region known as the critical layer (see Appendix B). However, we emphasise that, if one does not need precise order comparisons, reading § 2.3 is sufficient to follow the derivation of the PCS approach.

Section 2.5 outlines the relationship between the PCS and various asymptotic and parabolised approaches. Readers uninterested in these details may also skip this section.

2.2. Derivation of BRE

We take the Reynolds number $R$ as the large parameter in the asymptotic analysis. The rescaled variables suitable for the BRE are defined as $X=R^{-1}x$ and $T=R^{-1}t$ . The effect of varying curvature can be incorporated by treating $De$ as a function of $X$ . Substituting the asymptotic expansion

(2.3a) \begin{align} [u,v,w] &=[U,R^{-1}V,R^{-1}W](X,y,z,T)+\boldsymbol{\cdots} , \end{align}

(2.3b) \begin{align} p &=\varPi (X,T) +R^{-2}P(X,y,z,T)+\boldsymbol{\cdots} , \end{align}

into the governing equations, (2.1), and retaining only the leading-order terms yields the BRE

(2.4a) \begin{eqnarray} \kern-3.4pc \big[\partial _T+U\partial _X +V\partial _y +W\partial _z-\partial _y^2-\partial _z^2 \big] \left [\! \begin{array}{c} U \\ V\\ W \end{array} \!\right ] + \left [\! \begin{array}{c} \partial _X\varPi \\ \partial _y P+De^2U^2/8 \\ \partial _z P \end{array} \!\right ] &=0, \nonumber \\ \kern-3.4pc \partial _XU+\partial _yV+\partial _zW &=0. \end{eqnarray}

The coefficients $U,V,W,\varPi ,P$ appearing in the asymptotic expansion (2.3) are assumed to be $O(R^0)=O(1)$ , and the dots represent higher-order terms. Note that the substitution of the expansions and the subsequent order-of-magnitude analysis must be performed using the body-fitted coordinate form of the governing equations (see (A1) in the Appendix A). These algebraic manipulations are rather cumbersome, so rather than comparing the order of every term, we highlight only the key aspects.

The BRE system (2.4) closely resembles the Navier–Stokes equations for channel flow in Cartesian coordinates. Recalling that the walls are located at $y=\pm 1$ , the differences from the plane Poiseuille flow problem are: the presence of the term proportional to $De^2$ , the absence of the second derivative in $X$ and the absence of the $X$ derivative of $P$ . Owing to the latter two properties, the system is parabolic, and marching in $X$ becomes feasible.

Here, we briefly illustrate why the above features of the BRE arise. For clarity, we restrict the discussion to the constant-curvature case, where cylindrical coordinates $(r,\varphi ,z)$ can be used, since body-fitted coordinates are less familiar to most readers (see Appendix A for the general case). The well-known Navier–Stokes equations in cylindrical coordinates, expressed in terms of $\boldsymbol{u}=u_r\boldsymbol{e}_r+u_{\varphi }\boldsymbol{e}_{\varphi }+u_z\boldsymbol{e}_z$ , can be related to the body-fitted coordinate version by setting $u_r=-v, u_{\varphi }=u, u_z=w$ , $y=a-r$ and $x=a\varphi$ . The channel walls are located at $r=a\pm 1$ , with $a=O(R^2)$ , as seen in (2.2). This implies that $r=O(R^2)$ because within the channel $y\in (-1,1)$ .

We first note that $r^{-1}\partial _{\varphi }=r^{-1}aR^{-1}\partial _X\approx R^{-1}\partial _X$ . Thus the advection operator $u_r\partial _r+r^{-1}u_{\varphi }\partial _{\varphi }+u_z\partial _z \approx R^{-1}(U\partial _X +V\partial _y +W\partial _z$ ) is formally of size $O(R^{-1})$ , and may balance with the most terms in the viscous operator. Indeed, terms such as $R^{-1}(\partial _r^2+\partial _z^2+\ldots )\approx R^{-1}(\partial _y^2+\partial _z^2)$ have effects that are of the same order of magnitude, although second derivatives with respect to $X$ are negligible as $R^{-1}r^{-2}\partial _{\varphi }^2\approx R^{-3}\partial _X^2$ . The velocity scaling in (2.3a ) indicates that, in the momentum equations, the leading-order terms in the azimuthal component are $O(R^{-1})$ , whereas those in the radial and axial components are $O(R^{-2})$ . The pressure term in the azimuthal component can be approximated as $r^{-1}\partial _{\varphi }p\approx R^{-1}\partial _X\varPi +R^{-3}\partial _XP$ , and the second term is negligible. Since $r$ is large, most of the curvature-related terms can be neglected. However, in the radial component, the term $r^{-1}u_{\varphi }^2\approx a^{-1}u^2\approx R^{-2}De^2U^2/8$ produces the leading-order effect. The argument above applies in a very similar way to the case of non-constant wall curvature.

2.3. Coupling BRE with fast-scale local ECS: the PCS equations

Now, let us consider incorporating the effects of the fast-scale coherent structures into the BRE. Borrowing ideas from multiple-scale asymptotic analysis (see Howard & Kopell Reference Howard and Kopell1977, for example), we introduce a fast-scale variable $\theta =R\varTheta (X,T)$ , and decompose the velocity and pressure fields as

(2.5a) \begin{align} \boldsymbol{u} &=\overline {\boldsymbol{u}}(X,y,z,T)+\tilde {\boldsymbol{u}}(X,\theta ,y,z,T), \end{align}

(2.5b) \begin{align} p &=\overline {p}(X,y,z,T)+\tilde {p}(X,\theta ,y,z,T). \end{align}

The tilde components are assumed to be $2\pi$ periodic in $\theta$ . In what follows, the overline also denotes averaging over $\theta$ ; thus the above ‘Reynolds decomposition’ represents a separation into the $\theta$ -mean field and the fluctuation field. We will also decompose the governing equations into mean and fluctuation parts with respect to $\theta$ .

The assumed periodicity in $\theta$ is indeed consistent with the multiple-scale asymptotic analysis in § 2.4, but we emphasise that no asymptotic expansion is used in this section. Nevertheless, by formally considering large $R$ with $\varTheta = O(1)$ , it becomes evident that $\theta$ represents the fast scale, as shown by the following relation:

(2.6) \begin{eqnarray} \frac {\partial }{\partial t}=-\varOmega \frac {\partial }{\partial \theta } + R^{-1}\frac {\partial }{\partial T},\qquad \frac {\partial }{\partial x}=\alpha \frac {\partial }{\partial \theta } + R^{-1}\frac {\partial }{\partial X}, \end{eqnarray}

where we have written $\varOmega =-\varTheta _T$ and $\alpha =\varTheta _X$ . At this stage, the ‘local’ streamwise wavenumber $\alpha$ and frequency $\varOmega$ of the wave-like motion on the fast scale may depend on both $X$ and $T$ , provided they satisfy the consistency condition $\alpha _T+\varOmega _X=0$ . If $\alpha$ is taken as a function of $X$ and $\varOmega$ as a function of $T$ , then the consistency condition is automatically satisfied, and we have

(2.7) \begin{eqnarray} \varTheta (X,T)=\int ^X_0 \alpha (s){\textrm{d}}s-\int ^T_0 \varOmega (\tau ) {\textrm{d}}\tau . \end{eqnarray}

The mean field is independent of the fast scale and therefore physically represents fluid motion that evolves slowly on the BRE spatio-temporal scales.

To reduce the governing equations and facilitate numerical computation, we introduce the following assumptions.

1. (i) Assumption (i): the governing equations, (2.1), (equivalently, (A1)) can be simplified to

   (2.8a) \begin{align} \big(\partial _t +u\partial _x+v\partial _y+w\partial _z-R^{-1}\big(\partial _x^2+\partial _y^2+\partial _z^2\big)\big) \left [\! \begin{array}{c} u\\ v\\ w \end{array} \!\right ]& +\left [\! \begin{array}{c} \partial _x p\\ \partial _y p+u^2/a\\ \partial _z p \end{array} \!\right ]=0, \end{align}
   (2.8b) \begin{align} \partial _xu+\partial _yv &+\partial _zw=0. \end{align}

   We substitute the decomposition (2.5) into the governing equations and then separate them into mean and fluctuation parts by taking the $\theta$ -average.

2. (ii) Assumption (ii): the mean equations can be parabolised by neglecting $\partial _X$ in the viscous and pressure terms.

3. (iii) Assumption (iii): terms involving slow-scale derivatives and curvature effects in the Reynolds stress terms can be omitted.

4. (iv) Assumption (iv): likewise, terms involving slow-scale derivatives and curvature effects in the fluctuation equations can also be neglected.

We briefly explain here why they are expected to hold (the order of magnitude of the flow field is presented in § 2.4, while a more direct justification of the assumptions is given in Appendix A). Assumption (i) implies that, when the effect of curvature is weak, the flow may be approximated by the Cartesian Navier–Stokes equations, although the curvature terms appearing in the BRE must still be retained. Assumption (ii) has been confirmed for the BRE (§ 2.2). Assumptions (iii) and (iv) rely on the expectation that, when derivatives act on the fluctuation field, the contributions from the slow derivatives can be neglected; see (2.6).

Under the assumptions (i)–(iv), the $\theta$ -mean part of the governing equations reduce to the form of the forced BRE upon writing $\overline {\boldsymbol{u}}=[U,R^{-1}V,R^{-1}W]$ , $\overline {p}=\varPi (X,T) +R^{-2}P(X,y,z,T)$

(2.9a) \begin{align} \big[\partial _T+U\partial _X +V\partial _y +W\partial _z-\partial _y^2-\partial _z^2\big] \left [\! \begin{array}{c} U\\ V\\ W \end{array} \!\right ] + \left [\! \begin{array}{c} \partial _X\varPi \\ \partial _y P+De^2U^2/8 \\ \partial _z P \end{array} \!\right ] &=\boldsymbol{F}\!, \end{align}

(2.9b) \begin{align} \partial _XU+\partial _yV+\partial _zW &=0. \end{align}

This is not a surprising result, as in § 3.1 we have checked the first two assumptions are valid for spatio-temporally slowly developing flow. When assumption (iii) holds, the forcing term can be compactly written as

(2.9c) \begin{equation} \boldsymbol{F}=- \big[R\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\tilde {u}},R^2\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\tilde {v}},R^2\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\tilde {w}}\,\big], \end{equation}

which can be interpreted as the subgrid-scale Reynolds stress; here, $\tilde {\boldsymbol{\nabla }}=[\alpha \partial _{\theta },\partial _y,\partial _z]$ . As seen in § 2.2, obtaining the BRE requires rescaling the governing equations, and as a consequence several factors of $R$ appear in (2.9c ).

While the fluctuation part of the governing equations is obtained by subtracting the mean part from the full equations

(2.10a) \begin{align} (\overline {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }}-c\alpha \partial _{\theta })\tilde {\boldsymbol{u}}+ (\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\overline {\boldsymbol{u}} +(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}-\overline {(\tilde {\boldsymbol{u}}\boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\tilde {\boldsymbol{u}}} &=-\tilde {\boldsymbol{\nabla }} \tilde {p}+R^{-1}\tilde {{\nabla} }^2 \tilde {\boldsymbol{u}}, \end{align}

(2.10b) \begin{align} \tilde {\boldsymbol{\nabla }}\boldsymbol{\cdot }\tilde {\boldsymbol{u}} &=0. \end{align}

Here, $c(X,T)=\varOmega (T)/\alpha (X)$ is the local phase speed. Thanks to assumption (iv), no derivatives with respect to $X$ and $T$ appear in the fluctuation equations. Equations (2.9) and (2.10) constitute the PCS system we numerically solve to obtain the mean and fluctuation fields appearing in the decomposition (2.5). Note that the rescaling of $[\overline {\boldsymbol{u}},\overline {p}]$ with $U,V,W,\varPi ,P$ is merely formal in this section and is intended to emphasise that (2.9) corresponds to the forced BRE.

When $De = 0$ , the travelling-wave solutions of plane Poiseuille flow with phase speed $c$ exactly satisfy (2.9)–(2.10). This follows from the fact that $\alpha$ becomes constant and we can use the fast-scale variable $\theta =\alpha (x-ct)$ . The travelling waves fit within the periodic box $[0,2\pi /\alpha ] \times [-1,1] \times [0,2\pi /\beta ]$ in the $(x,y,z)$ coordinates, and therefore the dependence on the slow variables $X$ and $T$ disappears from the global flow. Even if $De$ is a non-zero constant, travelling-wave-type ECS solutions persist in the reduced system, although only for the approximate system (2.8a ) and not for the full cylindrical coordinate system.

The PCS (2.9)–(2.10) can thus be viewed as a means of incorporating Reynolds stresses generated by the SSP-VWI into the BRE via ECS. The asymptotic theory in § 2.4 and Appendix B will reveal that the fluctuation velocity $\tilde {\boldsymbol{u}}$ is slightly larger or smaller than $O(R^{-1})$ depending on the location, and it clarifies how the Reynolds stress terms balance with the BRE. Alternatively, the formulation may be regarded as a generalisation of ECS studies, which have so far been restricted to parallel flows, to include the spatial evolution of the system. Importantly, this is not merely a local parallel approximation; rather, it naturally accounts for what are referred to as non-parallel effects in boundary-layer theory.

2.4. Large-Reynolds-number asymptotic analysis

The assumptions (i)–(iv) used to derive the PCS equations can be justified using a large-Reynolds-number asymptotic analysis. Specifically, our argument here builds on the vortex/Rayleigh-wave interaction theory developed by Hall & Smith (Reference Hall and Smith1991), in which a multiple-scale asymptotic analysis is applied to boundary-layer flows. As noted in § 1, Hall & Sherwin (Reference Hall and Sherwin2010), Deguchi & Hall (Reference Deguchi and Hall2014b ), Gibson & Brand (Reference Gibson and Brand2014) and Deguchi & Hall (Reference Deguchi and Hall2016) confirmed that ECS in parallel flows follow the scaling inferred by this asymptotic theory.

In the asymptotic analysis, after the decomposition (2.5), the mean component is assumed to follow an asymptotic expansion similar to that in the BRE

(2.11a) \begin{align} \overline {\boldsymbol{u}} &=[U,R^{-1}V,R^{-1}W](X,y,z,T)+\boldsymbol{\cdots} , \end{align}

(2.11b) \begin{align} \overline {p} &=\varPi (X,T) +R^{-2}P(X,y,z,T)+\boldsymbol{\cdots} , \end{align}

while the fluctuation component is expanded as

(2.12a) \begin{align} \tilde {\boldsymbol{u}} &=\delta ^{1/2} R^{-1}e^{i\theta }\hat {\boldsymbol{u}}(X,y,z,T)+\text{c.c.}+\boldsymbol{\cdots} , \end{align}

(2.12b) \begin{align} \tilde {p} &=\delta ^{1/2} R^{-1}e^{i\theta }\hat {p}(X,y,z,T)+\text{c.c.}+\boldsymbol{\cdots} , \end{align}

where $\delta =R^{-1/3}$ and c.c. stands for complex conjugate. In the asymptotic expansions, the coefficients $U,V,W,\varPi ,P,\hat {\boldsymbol{u}},\hat {p}$ are all $O(1)$ . Consequently, the streamwise and cross-streamwise components of $\overline {\boldsymbol{u}}$ , called the streak ( $\bar u$ ) and roll ( $\bar v, \bar w$ ) in SSP-VWI, scale as $O(R^0)$ and $O(R^{-1})$ , respectively. The fluctuation part $\tilde {\boldsymbol{u}}$ , referred to as the wave, scales as $O(R^{-7/6})$ .

Substituting the asymptotic expansions (2.11)–(2.12) into (2.1) and retaining only the leading-order terms yields the BRE (2.4) together with the following equations:

(2.13a) \begin{align} i\alpha (U-c)\hat {u}+\hat {v}\partial _y U+\hat {w}\partial _z U &=-i\alpha \hat {p}, \end{align}

(2.13b) \begin{align} i\alpha (U-c)\hat {v} &=-\partial _y \hat {p}, \end{align}

(2.13c) \begin{align} i\alpha (U-c)\hat {w} &=-\partial _z \hat {p}, \end{align}

(2.13d) \begin{align} i\alpha \hat {u}+\partial _y\hat {v}+\partial _z\hat {w} &=0. \end{align}

These equations can be combined into a single equation for $\hat {p}$

(2.14) \begin{eqnarray} \partial _y\left (\frac {\partial _y\hat {p}}{(U-c)^2} \right )+\partial _z\left (\frac {\partial _z\hat {p}}{(U-c)^2} \right )-\alpha ^2 \frac {\hat {p}}{(U-c)^2}=0, \end{eqnarray}

namely the Rayleigh equation, which becomes singular when $U=c$ . This implies that the asymptotic expansions (2.11)–(2.12) break down near the critical level $y=f(z,X,T)$ , where $U=c$ . More specifically, this breakdown occurs within a thin layer of thickness $\delta =R^{-1/3}$ surrounding the critical level, known as the critical layer. Within this layer, a new asymptotic expansion must be applied, using an appropriate stretched coordinate. After a lengthy and involved algebra given in Appendix B, one concludes that the following conditions must be imposed on the mean component:

(2.15) \begin{eqnarray} \,[\partial _n (\boldsymbol{e}_s\boldsymbol{\cdot }(V\boldsymbol{e}_y+W\boldsymbol{e}_z))]^{f+}_{y=f-}=J(s,X,T),\qquad [P]^{f+}_{y=f-}=K(s,X,T). \end{eqnarray}

The left-hand sides represent the jumps across the critical level, while the right-hand sides are known functions that can be computed from $\hat {p}$ and $U$ (therefore, the mean part is influenced by the fluctuation part). Here, in the $y$ – $z$ plane, we introduce coordinates $n$ and $s$ , normal and tangent to the critical level, respectively ( $\boldsymbol{e}_s$ is the unit vector pointing in the $s$ direction). The leading-order system in the asymptotic theory is given by (2.4), (2.13) and (2.15). This is a closed system, but solving it is not easy due to the singularity at the critical level (see Deguchi Reference Deguchi2019, for example). Moreover, the location of that level is not known a priori as $U$ is unknown.

The important observation here is that substituting the asymptotic expansions inside and outside the critical layer into the PCS system (rather than into (A1)) yields the same leading-order system (2.4), (2.13) and (2.15). This clearly demonstrates that the PCS system retains all the leading-order terms in the asymptotic analysis. Appendix A provides a direct check of assumptions (i)–(iv), but only for selected terms from the governing equations, as a full verification would require considerable space.

We also provide a short summary of the main points of the analysis below. First, the equations for the wave component (2.13) has the form of the Navier–Stokes equations linearised about the streak $U$ . The reason for this is that, as already noted, the wave amplitude is smaller than $O(1)$ , while the streak component is $O(1)$ . In expansion (2.12), we assume that only a single Fourier mode in $\theta$ becomes neutral, as is the case for most known ECS, which explains why the leading-order terms are proportional to $e^{i\theta }$ . While the theory could, in principle, be extended to incorporate multiple neutral modes, such complications are not considered here. Those asymptotic structures of the wave can be obtained from the fluctuation equations in the PCS system (2.10), which retain higher-order wave–wave interaction terms.

Second, near the critical layer, viscous effects become significant, and the apparent singularity of the outer solution is thereby regularised. Outside the critical layer, advection $i\alpha (U-c)$ is $O(1)$ and dominates over viscous effects (e.g. $R^{-1}\partial _y^2=O(R^{-1})$ ); thus (2.13) is inviscid. Near the critical layer, however, the former becomes small, whereas the latter becomes significant on the short scale. The critical layer thickness $\delta \sim O(R^{-1/3})$ is determined by the balance of these two effects (see, for example, Drazin & Reid Reference Drazin and Reid2004). In PCS, by retaining viscosity, the structure inside and outside the critical layer can be described by a single set of equations.

Third, the jump condition (2.15) effectively accounts for the Reynolds stress within the critical layer. Due to the singularity, the fluctuation field $\tilde {\boldsymbol{u}}$ is amplified by a factor of $\delta ^{-1}$ there, yielding a magnitude of $O(\delta ^{-1/2}R^{-1})=O(R^{-5/6})$ . Thus, the size of the cross-streamwise components of the Reynolds stress $\boldsymbol{F}$ is amplified within the critical layer, balancing with the viscous terms in (2.9a ). Evaluating the effect of this mean–fluctuation interaction on the outer region leads to the jump conditions (2.15). Not surprisingly, the forcing term $\boldsymbol{F}$ in the PCS produces an equivalent effect to that of the jump conditions.

Finally, we note that the asymptotic expansion relies on a delicate balance, where even minor changes can lead to inconsistencies. For example, if a smaller wave size is chosen in (2.12), then $J=K=0$ in (2.15). In this case, the wave no longer feeds back to the roll–streak system, and the SSP-VWI cycle halts. Also, if the wave component is not neutral, the Reynolds stress grows or decays exponentially on the short scale, breaking the balance assumed in the asymptotic expansion. The choice of the phase function is also not completely arbitrary, as noted just below (2.6).

2.5. Additional remarks on the properties of the PCS system

The PCS system (2.9)–(2.10) retains more terms than the vortex–wave interaction theory. One advantage of the PCS approach is that it allows an ECS to be directly used as the inlet condition. Furthermore, when aiming to construct a good approximation of the Navier–Stokes solutions, it is desirable to minimise the number of discarded terms.

Another benefit is that the PCS can accommodate asymptotic solutions involving multiple Fourier modes. An example of such a case appears in the strongly nonlinear critical layer theory by Smith & Bodonyi (Reference Smith and Bodonyi1982) (for validation using ECS, see Deguchi & Walton Reference Deguchi and Walton2018). Also, when coherent structures appear in the free stream of the boundary layers (Deguchi & Hall Reference Deguchi and Hall2014a ) or at the Kolmogorov scale (Deguchi Reference Deguchi2015), all terms in the Navier–Stokes equations become important. Even in such extreme cases, assumptions (i)–(iv) remain valid, and the PCS equations can still be applied.

Hall & Smith (Reference Hall and Smith1991) also considered vortex/Tollmien–Schlichting wave interaction (see also Hall & Smith Reference Hall and Smith1988). A validation of this theory using ECS was later carried out by Dempsey et al. (Reference Dempsey, Deguchi, Hall and Walton2016) using plane Poiseuille flow. The PCS is capable of reproducing this type of asymptotic solution as well.

We also wish to comment on the differences between our method and another major spatial marching approach, the parabolised stability equations (PSE, see the seminal review by Herbert Reference Herbert1997). The PSE have been widely used to efficiently identify transition points to turbulence in various aerodynamic design problems (Chang Reference Chang2004). Rather than using (2.10) for the fast scale, PSE adopt a different equation that is not compatible with using an ECS as the inlet condition. A more fundamental difference lies in the treatment of the streamwise wavenumber $\alpha$ , which in PSE is allowed to be complex and is updated at each step of the spatial marching, with several standard iterative adjustment strategies available (see Bertolotti, Herbert & Spalart Reference Bertolotti, Herbert and Spalart1992). The imaginary part of $\alpha$ represents the spatial growth rate of the fluctuation field. Consequently, if this growth rate is non-zero, we cannot mathematically exclude the possibility that the Reynolds stresses vary on the fast scale. This violates a key assumption of multiple-scale analysis, namely that the mean field depends only on the slow spatial variables. Since the Reynolds stress is a nonlinear quantity, this issue arises when the perturbations have finite amplitude.

Although reasonably good agreement between PSE results and asymptotic theory has often been observed, this agreement breaks down over long spatio-temporal scales. Readers interested in asymptotic analyses of spatially developing flows may consult Wu (Reference Wu2019). It should also be noted that the objectives of PSE and our approach differ fundamentally. While the former seek to capture the evolution of the flow field near the transition point, the latter aims to uncover the key ECS that organise the phase space.

3.1. Computational set-up

In § 4, we numerically investigate a curved channel with the Dean number $De(X)=De_{\textit{max}}\,\chi (X)$ , where

(3.1) \begin{equation} \chi (X) = \begin{cases} 0, & 0 \leqslant X \lt X_{\!{A}}, \\ \left[\frac {1}{2}-\frac {1}{2}\cos \left(\frac {X-X_{\!{A}}}{X_{\!{B}}-X_{\!{A}}}\pi \right)\right]^{{1}/{2}} \!, & X_{\!{A}} \leqslant X \lt X_{\!{B}}, \\ 1,& X_{\!{B}} \leqslant X. \end{cases} \end{equation}

Our choices of the parameters are $De_{\textit{max}}=\sqrt {60\,500}\approx 245.9$ , $X_{\!{A}}=10$ and $X_{\!{B}}=130$ . The profile $De(X)$ is shown in figure 2(a). According to the stability analysis of Dean flow in a uniformly curved channel, an instability arises when $De$ exceeds 35.9 (see Gibson & Cook Reference Gibson and Cook1974; Mizushima, Yanase & Yoshizawa Reference Mizushima, Yanase and Yoshizawa1998, for example). Thus, $De_{\textit{max}}$ used in our study is large enough to induce Dean vortices. The flow in the region $X \lt 0$ is assumed to be plane Poiseuille flow, within which an edge state is sustained. We choose $R = 550$ and carry out spatial marching of the reduced (2.9)–(2.10) over the interval $X \in [0, 150]$ .

Figure 2. Channel curvature. (a) Dean number profile defined in § 3.1, used for spatial marching in §§ 4.1 and 4.2. (b) Channel centreline plotted in Cartesian coordinates.

The Dean number determines the radius of curvature of the channel centreline via (2.2). As shown in figure 2(b) the channel forms a spiral. Here, to visualise the actual channel geometry, we introduce Cartesian coordinates $(\check {x}, \check {y}, \check {z})$ , which coincide with the curvilinear coordinates $(x, y, z)$ at the origin; see figure 1 also. A total arc length of the spiral is $550\times 150=82\,500$ . Computing ECS over such a large domain using the full Navier–Stokes equations would be impractical. Towards the end of the spiral, the radius of curvature reduces to 40, which nevertheless remains much larger than the channel width.

Since $De(X)$ vanishes at the inlet $X=0$ , a travelling-wave-type ECS of plane Poiseuille flow can naturally serve as the initial condition for spatial marching. Note that such a travelling wave is independent of the slow time variable $T$ , and therefore the global solution of the reduced system also becomes independent of $T$ . The fast-scale variable simplifies

(3.2) \begin{eqnarray} \theta =R\int ^X_0 \alpha (s){\textrm{d}}s- \varOmega t, \end{eqnarray}

where $\varOmega$ is a constant.

In the following sections, $\overline {\boldsymbol{u}}$ depends only on $X,y,z$ , and is therefore steady. On the other hand, $\tilde {\boldsymbol{u}}$ also depends on $\theta$ defined in (3.2). The solution obtained by spatial marching is time periodic with period $2\pi /\varOmega$ ; recall $\theta$ is normalised so that the solution has a $2\pi$ period with respect to this variable.

The marching problem depends on the frequency $\varOmega$ and the spanwise wavenumber $\beta$ ; these must be determined so that they are consistent with the inlet ECS. The pressure gradient $\varPi '(X)$ and the local wavenumber $\alpha (X)$ are determined as part of the solution.

3.2. Inlet conditions: ECS of plane Poiseuille flow

We employ two types of travelling-wave solutions of plane Poiseuille flow, both obtained by solving the full Navier–Stokes equations in a periodic box $(x,y,z)\in [0,2\pi /\alpha ]\times [-1,1]\times [0,2\pi /\beta ]$ using the Newton–Raphson method. Both solutions are found via symmetry-breaking bifurcations from the MS-S solution reported by Nagata & Deguchi (Reference Nagata and Deguchi2013). To describe the symmetry of the solution, we use the fact that a travelling wave with phase speed $c$ is a function of $\theta$ , $y$ and $z$ , where

(3.3) \begin{eqnarray} \theta =\alpha (x-ct). \end{eqnarray}

The MS-S solution has the three symmetries : the mirror symmetry (MS) in the spanwise direction

(3.4) \begin{eqnarray} \,[u,v,w](\theta ,y,z)=[u,v,-w](\theta ,y,-z+\pi /\beta ), \end{eqnarray}

the shift–reflection symmetry

(3.5) \begin{eqnarray} \,[u,v,w](\theta ,y,z)=[u,v,-w](\theta +\pi ,y,-z), \end{eqnarray}

and top–bottom reflectional symmetry

(3.6) \begin{eqnarray} \,[u,v,w](\theta ,y,z)=[u,-v,w](\theta ,-y,z). \end{eqnarray}

Figure 3. Symmetry-breaking bifurcation from the MS-S solution branch in plane Poiseuille flow. Panels show (a) $(\alpha ,\beta )=(1.0,2.0)$ ; (b) $(\alpha ,\beta )=(1.01,2.72)$ .

Figure 4. Bifurcation diagram for plane Poiseuille flow, showing the pressure-based Reynolds number $R_p$ versus the bulk Reynolds number $R$ . The wavenumbers are optimised to minimise the value of $R_p$ for each solution. The solid lines represent the ECS. The solution names are provided in the legend, along with the corresponding wavenumbers $(\alpha , \beta )$ . The red and blue circles indicate the initial solutions used in the spatial marching shown in figures 7 and 9, respectively (see table 1 also).

Figures 3(a) and 3(b) illustrate the bifurcation diagram for two sets of wavenumber pairs $(\alpha ,\beta )$ . The vertical axis is $R_p=U_c^{*}d^{*}/\nu$ defined using the centreline velocity of laminar flow, $U_c^{*}$ ; this Reynolds number is proportional to the pressure gradient required to maintain the bulk velocity at unity. At the bifurcation point indicated by the green circle, the symmetry of MS-S in the wall-normal direction (3.6) is broken, while the other two symmetries are preserved. We refer to the solutions shown in panels (a) and (b) as MS-A and MS-A2, respectively (thus S and A represent top–bottom symmetry and asymmetry, respectively).

In figure 4 the wavenumber pairs $(\alpha ,\beta )$ are selected to minimise the onset $R_p$ of each solution. This bifurcation diagram, presented in nearly the same format as figure 5 of Nagata & Deguchi (Reference Nagata and Deguchi2013), clearly demonstrates that MS-A and MS-A2 are distinct solutions, despite sharing the same flow symmetry. As the Reynolds number increases, two solution branches emerge abruptly via a saddle node bifurcation. For the inlet conditions of the spatial marching problem presented in § 4, we use the lower-branch solutions of MS-A and MS-A2 at $R = 550$ (indicated by circles). The parameters of these solutions are summarised in table 1, and the corresponding flow fields are shown in figures 5(b) and 5(c). Note that the frequency $\varOmega$ of the global problem must be equal to the product $\alpha c$ determined from the inlet ECS. Also, the pressure gradient at $X=0$ is given by $\varPi '=-2R_p/R$ .

Table 1. Parameters and key physical properties of the initial plane Poiseuille ECS at the inlet, used in the spatial marching problem discussed in § 4. These solutions correspond to the circles in figure 4.

Figure 5. Flow visualisation of the lower-branch ECS in plane Poiseuille flow at $R=550$ . The optimal wavenumbers used in figure 4 are adopted here. The red/blue surfaces depict 50 % maximum/minimum value of $\partial _y\tilde {w} - \partial _z\tilde {v}$ . The brown/green surfaces are 50 % maximum/minimum value of $\partial _y\overline {w}-\partial _z \overline {v}$ . The colour map illustrates $\overline {u}$ . The thick red lines indicate the location of the critical levels, where $\overline {u}=c$ . In the absence of slow-scale dependence, $(\overline {u},\overline {v},\overline {w})$ correspond to the streamwise-averaged velocity fields.

It is known that the lower-branch solutions are well described by the vortex/Rayleigh-wave interaction theory, even at moderately high Reynolds numbers (Hall & Sherwin Reference Hall and Sherwin2010; Deguchi & Hall Reference Deguchi and Hall2014b ; Gibson & Brand Reference Gibson and Brand2014). Indeed, the flows shown in figure 5 exhibit all three elements of the SSP-VWI: rolls, streaks and waves. The presence of waves and rolls can be identified from the isosurfaces of streamwise vorticity in the fluctuation and mean fields, respectively. The streak field can be seen in the colour map. The two red curves indicate the critical levels, around which a strong wave field is generated, as predicted by Hall & Sherwin (Reference Hall and Sherwin2010).

Another important feature of the lower-branch solution is that, when the Navier–Stokes equations are viewed as a dynamical system, it lies on the basin boundary separating the laminar and turbulent attractors in the phase space (Itano & Toh Reference Itano and Toh2001; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007; Wang et al. Reference Wang, Gibson and Waleffe2007; Deguchi & Hall Reference Deguchi and Hall2016). Our inlet conditions correspond to an edge state in the sense that small perturbations around them can lead either to turbulence or relaminarisation. (Note that some studies define edge states more narrowly, referring specifically to the converged solution obtained through the flow control algorithm known as edge tracking.)

3.3. Numerical scheme

To solve (2.9)–(2.10) numerically, we write the mean field and the fluctuation field as

(3.7) \begin{eqnarray} \left [\! \begin{array}{c} U\\ V\\ W \end{array} \!\right ]= \left [\! \begin{array}{c} -\varPsi _z+\varXi _y \\ -\varPhi _{zz}-\varXi _X\\ \varPhi _{yz}+ \varPsi _{X}\end{array} \!\right ]\!,\qquad \left [\! \begin{array}{c} \tilde {u}\\ \tilde {v}\\ \tilde {w} \end{array} \!\right ]= \left [\! \begin{array}{c} \alpha \phi _{\theta y}-\psi _z \\ -\alpha ^2\phi _{\theta \theta }-\phi _{zz}\\ \phi _{yz}+\alpha \psi _{\theta } \end{array} \!\right ]\!, \end{eqnarray}

where partial derivatives have been expressed using subscripts to simplify the notation. As a result, the primitive variables $[U, V, W, P, \tilde {u}, \tilde {v}, \tilde {w}, \tilde {p}]$ can be reduced to a smaller set of variables $[\varXi , \varPhi , \varPsi ,\phi , \psi ]$ .

Here, $\phi (X,\theta ,y,z)$ and $\psi (X,\theta ,y,z)$ represent the poloidal and toroidal potentials of the fluctuation field, respectively. The use of such potentials is well established in computations with periodic boxes (see Schmitt & Von Wahl Reference Schmitt, Von Wahl, Heywood, Masuda, Rautmann and Solonnikov1992 and McBain Reference McBain2005 for example), and it naturally extends to the fluctuation field in our formulation. The governing equations for $\phi$ and $\psi$ can be found by applying the operators $\boldsymbol{e}_y \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \times \tilde {\boldsymbol{\nabla }} \times$ and $\boldsymbol{e}_y \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \times$ to (2.10a ), respectively. The continuity is automatically satisfied in this formulation.

A similar potential decomposition can also be applied to the mean field. However, for the BRE, the $z$ -averaged components require special treatment, as discussed by Deguchi et al. (Reference Deguchi, Hall and Walton2013). Thus, we decompose the field into its $z$ -average and the remainder, and introduce the potentials $\varPhi (X,y,z)$ and $\varPsi (X,y,z)$ for the latter, ensuring they have zero $z$ -average. Then, for the components $\langle U \rangle$ and $\langle V \rangle$ , where angle brackets denote averaging in the $z$ direction, we introduce a streamfunction $\varXi (X,y)$ . If the inlet condition satisfies MS in the $z$ -direction, the global solution inherits this symmetry, and thus $\langle W \rangle =0$ . The pressure $P$ can be eliminated from the $z$ fluctuation part of (2.9a ), yielding two equations for $\varPhi$ and $\varPsi$ . The equation for $\varXi$ is obtained by taking the $z$ -average of the streamwise component of (2.9a )

(3.8) \begin{eqnarray} \partial _X\langle U^2\rangle +\partial _y\langle UV\rangle -\partial _y^2 \langle U\rangle +\varPi '=\langle \boldsymbol{e}_x\boldsymbol{\cdot }\boldsymbol{F}\rangle . \end{eqnarray}

The equations for $\Phi$ and $\Psi$ , together with (3.8), are equivalent to the forced BRE and can be integrated in $X$ by applying finite differences.

If the fluctuation field is known, the problem reduces to solving the BRE with extra forcing terms. Conversely, if the mean field is known, the equations for $\phi$ and $\psi$ can be solved by the Newton–Raphson method, in much the same way as in periodic box ECS computations. Thus, the most straightforward strategy for solving (2.9)–(2.10) is to iteratively update the mean and fluctuation fields at each spatial marching step until convergence is achieved. However, unfortunately, convergence is not guaranteed in general. Therefore, in this study, we solve both fields simultaneously to ensure robust convergence. Upon applying backward differences to the $X$ derivatives, we obtain implicit inhomogeneous equations for the updated fields. This system can be efficiently solved using the Newton–Raphson method, by extending the code developed by Deguchi et al. (Reference Deguchi, Hall and Walton2013).

The Fourier spectral method is used in the $\theta$ and $z$ directions, while the Chebyshev collocation method is employed in the $y$ direction. We expand $\varXi$ , $\varPhi$ , and $\phi$ in the basis $(1 - y^2)^2 T_l(y)$ , and $\varPsi$ and $\psi$ in the basis $(1 - y^2) T_l(y)$ so that the no-slip conditions

(3.9) \begin{eqnarray} \varXi =\partial _y\varXi =\varPhi =\partial _y\varPhi =\phi =\partial _y\phi =\varPsi =\psi =0\qquad \text{at} \qquad y=\pm 1 \end{eqnarray}

are fulfilled. Here, $T_l(y)$ is the $l$ th Chebyshev polynomial, and the collocation points are $y_k=\cos ((k+1)\pi /(K+2)), k=0,\ldots ,K$ .

The choice of basis functions for $\varXi$ imposes four conditions that ensure $\langle U\rangle =\langle V\rangle =0$ on both walls. Thus integrating continuity $\langle U\rangle _X+\langle V\rangle _y=0$ with respect to $y$ over $[-1,1]$ , we see that the flux $\int ^1_{-1}\langle U\rangle {\textrm{d}}y$ is unchanged in $X$ . To obtain a fourth-order differential equation for $\varXi$ , we use the equation obtained by differentiating (3.8) with respect to  $y$ . The eliminated $\varPi '$ can be recovered using $\varPi ' = (\partial _y^2 \langle U \rangle ) |_{y = -1}$ , which follows from evaluating (3.8) at the lower wall. It can also be computed by integrating (3.8) in $y$

(3.10) \begin{eqnarray} \partial _X\int ^1_{-1}\langle U^2\rangle {\textrm{d}}y-[\partial _y \langle U\rangle ]^1_{-1}+2\varPi '=0. \end{eqnarray}

We analytically compute the Jacobian matrix and solve the linear system arising in each Newton step using a direct method. To eliminate the translational degree of freedom in the $\theta$ direction, we set the imaginary part of one spectral coefficient associated with the $e^{i\theta }$ mode to zero. The converged solution of the Newton iterations is obtained by treating $c$ as an unknown while keeping $\alpha$ fixed. To match $\alpha c$ with the given $\varOmega$ , a line search is used to find the optimal $c(\alpha )$ that minimises $|1-{\alpha c}/{\varOmega }|$ to $O(10^{-7})$ . The convergence of the search is fast, so each spatial marching step incurs a cost comparable to that of computing an ECS in a small periodic domain.

In the absence of the fluctuation field, the code is validated against the nonlinear Görtler vortex solution obtained by Hall (Reference Hall1988). We also verified that, when $De=0$ , the solution for plane Poiseuille flow presented in § 3.2 can be reproduced.

Figure 6. Key quantities obtained from the spatial marching problem in the curved channel. The red and blue curves represent the cases where MS-A and MS-A2 in table 1 are used as the inlet conditions, respectively. (a) Local streamwise wavenumber $\alpha$ . At $X=88.35$ , the blue curve terminates because the fast-scale dependence disappears and the Dean vortex solution emerges. (b) Pressure gradient $-\varPi '$ . The dot-dashed curve corresponds to the Dean vortex solution.

4.1. Marching with MS-A as the inlet condition

We begin by using MS-A from table 1 as the inlet condition and describe the typical outcome of the spatial marching. As noted in § 3.1, both $\alpha$ and $\varPi '$ are obtained as part of the solution. These results are shown as the red curves in figure 6.

The channel-curvature profile, given in (3.1), is designed to include a straight section over the range $X \in (0,10)$ . In this region, the curvature is zero, and the quantities $\alpha$ and $\varPi '$ , plotted in figure 6, remain constant. This confirms that MS-A is stable under spatial marching in the plane Poiseuille flow, despite being unstable in direct numerical simulations. As the curvature gradually increases along the section $X \in [10,130)$ , both $\alpha$ and $\varPi '$ undergo non-trivial changes. Once the curvature becomes constant again beyond $X = 130$ , these quantities rapidly converge to steady values.

Figure 7. The flow fields obtained for the curved channel flow using the spatial marching method. The channel curvature is defined by $De(X)$ with $De_{\textit{max}}=\sqrt {60\,500}\approx 245.9$ , as shown in figure 2. The inlet condition at $X=0$ is the lower branch of MS-A (see table 1). (a) Mean field. Isosurfaces show the roll component, and the colour map illustrates the streak component. (b) Snapshots of the fluctuation field at the locations corresponding to the colour map positions in panel (a). Panels (c) and (d) show the same plots as (b), taken at one quarter and one half of a time period later, respectively. See the caption of figure 5 for precise definition of the isosurfaces.

The numerical computation also determines the mean field, $\overline {\boldsymbol{u}}$ , and the fluctuation field, $\tilde {\boldsymbol{u}}$ , as defined in (2.5a ), at each streamwise location $X$ . Since their spatial scales differ significantly, we choose to plot them separately in figure 7. Plotting the mean (slow-scale) field in the $(X,y,z)$ coordinates is straightforward, as shown in panel (a). The isosurfaces represent the streamwise vorticity, $\partial _y\overline {w} - \partial _z\overline {v}$ , while the colour map shows the streamwise velocity, $\overline {u}$ . These quantities are steady and correspond to the roll and streak components of the SSP-VWI, respectively (note that our solutions do not depend on the slow time $T$ , as remarked in § 3.1). In contrast, plotting the fluctuation (fast-scale) field requires computing $\theta$ using (3.2). For instance, when fixing $t = 0$ , the value of $\theta$ at each $X$ is determined by integrating the streamwise wavenumber $\alpha$ shown in figure 6(a). Figure 7(b–d) shows the snapshots of the streamwise vorticity, $\partial _y\tilde {w} - \partial _z\tilde {v}$ ; they visualise the wave component in the SSP-VWI. Since the fluctuation field has a short characteristic streamwise length scale of $x=O(1)$ , we use the small boxes defined around $X=X_0=10, 40, 70, 100$ and $130$ . The boxes have dimensions $[0,2\pi /\alpha (X_0)] \times [-1,1] \times [0,2\pi /\beta ]$ in the shifted original coordinates $(x-RX_0,y,z)$ .

The fluctuation fields in panels (b–d) appear nearly as travelling waves, but not exactly. This is because both the fluctuation field functions $\tilde {\boldsymbol{u}}$ and the phase $\theta$ vary with the streamwise position. Note also that the boxes in panels (b–d) maintain the correct aspect ratio in the original coordinates $(x,y,z)$ . When converted to the $X$ -coordinate, the local wavelength of the fluctuation wave is only $O(0.01)$ .

Figure 7 shows that the centrifugal force induced by the channel curvature reinforces the mean field in the lower half of the channel, consistent with the findings of Dean (Reference Dean1928) and subsequent studies. However, there is no sharp boundary between the self-sustained coherent structures and the Dean vortices, as the fluctuation part of the solution persists throughout the spatial marching. At the inlet, wave-like structures associated with the fluctuation field are present near both walls, but those in the upper half of the channel are gradually suppressed. This suppression occurs because the weakening of the roll–streak structure by the centrifugal force disrupts the SSP-VWI cycle.

Figure 8. Comparison of constant-curvature channel flow solutions obtained from the reduced and full systems. (a) Solution obtained from the spatially marching computation shown in figure 7. This solution is taken near the outlet at $X=150$ (where $De=De_{\textit{max}}=245.9$ ) and is plotted in a local small box, similar to figure 7(b). See the caption of figure 5 for precise definitions of the isosurfaces and colour maps. (b) The corresponding travelling-wave solution of the Navier–Stokes equations in the annulus, with the bulk Reynolds number, frequency and radius of curvature matching those of (a) (i.e. $(R, \varOmega , a) = (550, 1.44, 40)$ ). For the relationship between the body-fitted and cylindrical coordinates, see § 2.2. Panels (c) and (d) show the same plots as (a) and (b), respectively, but for $De_{\textit{max}} = 50$ .

At the outlet, the slow-scale effects have essentially vanished, and the solution approaches a travelling wave in a constant-curvature channel. This is evidenced in figure 6, where the curves become nearly flat at $X=150$ . Since we are solving the reduced system (2.8a ), it is of interest to assess the accuracy of this approximation. To verify whether it provides a good approximation to the ECS of the full governing equations, (2.1), we employ the numerical code developed by Deguchi & Nagata (Reference Deguchi and Nagata2011), which computes ECS in concentric cylinders. The corresponding travelling-wave solution of the full Navier–Stokes equations can be obtained in cylindrical coordinates by considering flow between two concentric cylinders with a mean radius of curvature $a=40$ .

Figure 8(a) shows the flow field of the outlet solution from the reduced system, while panel (b) presents that of the ECS solution in the cylindrical Navier–Stokes equations. Comparison is made at the same bulk Reynolds number and frequency. Overall, the agreement is reasonable. The outlet solution has $\alpha =1.41$ , which compares well with $1.44$ for the full Navier–Stokes solution (in cylindrical coordinates, $\alpha$ is defined as the azimuthal wavenumber multiplied by $a$ ). Note that some differences in the flow details are not surprising, as $De_{\textit{max}}$ was chosen near the limit where $|\varPi '|$ remains close in value (less than 5 %) for both ECS. If $De_{\textit{max}}$ is reduced, the agreement improves significantly, as shown in panels (c) and (d), because the small-channel-curvature assumption is more easily satisfied. At this level of curvature ( $a=O(10^3)$ ), the errors in $(\alpha ,-\varPi ')$ are less than 0.1 %.

4.2. Marching with MS-A2 as the inlet condition

A qualitatively different result arises in the spatial marching problem when MS-A2 from table 1 is used as the inlet condition. The values of $\beta$ and $\varOmega$ are set to 2.82 and 1.52, respectively, as listed in the table, while all other computational conditions are kept the same as in § 4.1.

Figure 9. Same computation and plot as figure 7, but with the inlet initial condition set to MS-A2 (the blue circle in figure 4, also see table 1). Only a snapshot of the fluctuation field at one representative time is shown for simplicity. We use the endpoint value in figure 6(a), $\alpha = 1.51$ , to visualise the two rightmost boxes in panel (b).

Figure 10. Detailed view of the flow field near the location where the fluctuation field vanishes in figure 9. (a) Only the mean field is shown. The brown/green surfaces depict 50 % maximum/minimum of the streamwise vorticity, $\partial _y\overline {w}-\partial _z \overline {v}$ . (b) Enlarged view of the small box indicated in panel (a), with the fluctuation field vorticity overlaid (red/blue isosurfaces are 10 % maximum/minimum of $\partial _y\tilde {w}-\partial _z \tilde {v}$ ).

The result is shown in figure 9. Panel (a) shows the mean field, which consists of the roll and streak structure enhanced by centrifugal instability; this is qualitatively similar to the previous case. However, as seen in panel (b), the fluctuation field disappears before reaching $X=100$ . This transition occurs around $X = 88.35$ , as highlighted in the close-up view (figure 10). Beyond this critical point, the computation reduces to the BRE formulation (2.4), and the value of $\alpha$ does not affect the flow field (see figure 6 a).

The solution for $X \gt 88.35$ can be interpreted as a Dean vortex solution, which does not depend on the fast scale. The mean field at the outlet $X = 150$ satisfies the approximate (2.8a ) with $\partial _x=0$ and constant $De$ , which is none other than the reduced system originally derived by Dean. The comparison with the full Navier–Stokes solution in the annulus, shown in figure 11, confirms that the approximation is reasonably accurate, with the value of $|\varPi '|$ deviating by only 1 %.

The blue dot-dashed curve in figure 6(b) shows that emergence of Dean vortices increases the pressure gradient $|\varPi '|$ . Nevertheless, it remains consistently smaller than the red curve obtained using MS-A. The inlet–outlet pressure difference driving the latter flow is 578.98, whereas that of the former is significantly lower at 496.06.

Figure 11. Comparison between the reduced system solution and the Navier–Stokes solution in the constant-curvature case corresponding to the case shown in figure 9. The format of the figures is the same as figure 8. (a) Shows the Dean vortex solution taken near the outlet $X=150$ . (b) Shows the steady axisymmetric Navier–Stokes solution in the annulus, maintaining the same bulk Reynolds number.

Figure 12. The results of spatial marching computation with the inlet and outlet of the curved channel reversed. (a) The outcome of the spatial marching computation with the inlet condition chosen as the Dean vortex solution shown in figure 11(a). Along the dot-dashed line the Dean vortices are sustained, while the dashed line represents the unidirectional laminar solution. The green solid line corresponds to the result with a small external forcing of frequency 1.518. (b) Growth rate obtained from the secondary stability analysis of Dean flow. (c) The corresponding frequency of the eigenmode.

4.3. Transition from Dean vortices to self-sustained states

In § 4.2, the numerical computation was initiated from a self-sustained ECS of plane Poiseuille flow, MS-A2, and resulted in a Dean vortex solution at the outlet, $X = 150$ . A natural question then arises: What happens when the Dean vortex solution is used as the initial condition in a system with reversed inlet and outlet? Is it possible to recover MS-A2 after the channel becomes straight?

The answer to this question is no, if the spatial marching is initiated from the Dean vortex solution shown in figure 11(a) (i.e. the solution at $X=150$ in figure 6(b), depicted by the blue dashed curve). The result of the spatial marching is shown by the blue line in figure 12(a). Here and hereafter, we switch to the curvature $De(X) = De_{\textit{max}}\,\chi (150 - X)$ , while keeping the inlet position at $X=0$ . As the curvature decreases with $X$ , the Dean vortices (dot-dashed line) gradually weaken and vanish beyond $X\approx 70$ (the unidirectional laminar flow, indicated by the dashed line, has $\varPi '=-3$ ). The results here show that while spatial marching readily captures simple, highly symmetric solutions at downstream, it struggles to identify more complex, low-symmetry solutions without additional care. The purpose of this section is to outline a general method for obtaining such low-symmetry solutions via spatial marching. The problem set-up above is convenient because the existence of MS-A2 at the outlet is already known.

To reproduce the MS-A2 state from the inlet Dean vortices as demonstrated by the green line in figure 12(a), only a subtle adjustment to the previous computation is needed. We first note that, when the fluctuation field is infinitesimally small, the fluctuation (2.10) can be linearised with respect to that field. Writing $\tilde {\boldsymbol{u}}=e^{i\theta }\hat {\boldsymbol{u}}(X,y,z)$ and $\tilde {p}=e^{i\theta }\hat {p}(X,y,z)$ , the linearised equations become

(4.1a) \begin{align} {\mathscr{L}}\hat {u}+\hat {v}\partial _y \overline {u}+\hat {w}\partial _z \overline {u} &=-i\alpha \hat {p}, \end{align}

(4.1b) \begin{align} {\mathscr{L}}\hat {v}+\hat {v}\partial _y \overline {v}+\hat {w}\partial _z \overline {v} &=-\partial _y \hat {p}, \end{align}

(4.1c) \begin{align} {\mathscr{L}}\hat {w}+\hat {v}\partial _y \overline {w}+\hat {w}\partial _z \overline {w} &=-\partial _z \hat {p}, \end{align}

(4.1d) \begin{align} i\alpha \hat {u}+\partial _y\hat {v}+\partial _z\hat {w} &=0, \end{align}

where ${\mathscr{L}}=i\alpha (\overline {u}-c)+\overline {v}\partial _y+\overline {w}\partial _z-R^{-1}(\partial _y^2+\partial _z^2-\alpha ^2)$ . The above set of equations represents the temporal secondary stability problem under the parallel-flow approximation, a well-established approach in boundary-layer theory; see Hall & Horseman (Reference Hall and Horseman1991), Li & Malik (Reference Li and Malik1995) and Xu et al. (Reference Xu, Zhang and Wu2017), for example (note that Hall & Horseman Reference Hall and Horseman1991 used (2.13)). At each $X$ , the Dean vortex is substituted into $\overline {\boldsymbol{u}}$ in (4.1), and the equations are solved for the eigenvalue $c=c_r+ic_i$ . The result with $\alpha =1.51$ is shown in figures 12(b) and 12(c); the value of $\alpha$ is taken from the point where the blue curve in figure 6(a) terminates. The location where the growth rate becomes zero is at $X=61.64$ , and the corresponding frequency at that point is 1.518.

Figure 13. (a) Results of applying external forcing with various frequencies to the same marching problem as in figure 12(a). (b) The shaded region indicates the linearly unstable regime predicted by the stability analysis of the Dean vortex solution. The horizontal lines correspond to the solutions shown in panel (a).

Figure 14. The MS-A3 solution. (a) Bifurcation diagram shown in the same format as figure 4. (b) Flow visualisation of the lower-branch MS-A3 solution at $R=550$ (circle in panel (a)). See the caption of figure 5 for definitions of the isosurfaces and the colour map.

We then add a small random external forcing term with $\theta$ dependence to (2.10a ). This forced marching problem is initiated from the Dean vortex at $X=60$ , which is slightly upstream of the neutral point. The value of forcing frequency $\varOmega =1.518$ is chosen to match the one obtained at the neutral point. The result corresponds to the green curve in figure 12(a) mentioned earlier. The solution at $X=150$ converges to MS-A2, although the value of $\alpha$ differs slightly from that shown in table 1.

If the local parallel approximation were used by setting $\partial _X=0$ , the above computation could be interpreted as an imperfect bifurcation, as the forcing breaks the $\theta$ -invariance of the Dean vortices. However, since we solve the evolution equations along the streamwise direction, this interpretation is not exact. Instead, the problem is more closely related to the receptivity problem in aeronautics (Morkovin Reference Morkovin1969; Saric, Reed & Kerschen Reference Saric, Reed and Kerschen2002), which concerns the relationship between the input (external forcing from wall roughness or free-stream perturbations) and the output (downstream perturbations) in a linear system. Near the neutral point, our analysis can be interpreted as a nonlinear extension of this framework. Another important distinction from conventional boundary-layer transition is that, in our case, the base flow (Dean vortex) is initially unstable in the upstream region, and receptivity becomes most effective when the flow stabilises. To the best of our knowledge, such ‘subcritical’ receptivity has not been explored before.

Similar computations can be performed for various forcing frequencies, as shown in figure 13(a). Panel (b) shows the frequency range where instability occurs. When the frequency $\varOmega$ falls below approximately 1.2, the marching terminates before reaching the outlet due to the failure of the Newton method to converge. Interestingly, however, once $\varOmega$ is decreased further to around 0.55, the marching suddenly becomes stable and can proceed all the way to the outlet.

The plane Poiseuille solution found at this frequency is qualitatively different from any of the solutions introduced in § 3.2, and we therefore refer to it as MS-A3. This difference can be seen from the bifurcation diagram in figure 14(a) and the flow field in figure 14(b), both presented in a format that allows direct comparison with figures 4 and 5, respectively. These computations demonstrate that the marching method can serve as a practical tool for discovering new ECS in parallel flows.