1. Introduction
In isothermal pipe flow, the flow is driven by an external pressure gradient. This is referred to as ‘forced’ flow. In a vertical configuration, however, buoyancy resulting from the expansivity of the fluid close to a heated wall can provide a force that partially or fully drives the flow, referred to as ‘mixed’ or ‘natural’ convection, respectively. The use of mixed convection is fundamentally important and practical in engineering applications, e.g. geothermal energy capture, nuclear reactor cooling systems, fossil fuel power plants, and has been widely researched. Heat transfer in the presence of buoyancy exhibits totally different characteristics for downward flow and upward flow. In a downward flow, the buoyancy acts as a drag force but always enhances the heat transfer. By contrast, in an upward flow, buoyancy assists the flow, but heat transfer can deteriorate significantly (Ackerman Reference Ackerman1970; McEligot, Coon & Perkins Reference McEligot, Coon and Perkins1970; Jackson & Li Reference Jackson and Li2002; Bae, Yoo & Choi Reference Bae, Yoo and Choi2005; Wibisono, Addad & Lee Reference Wibisono, Addad and Lee2015; Zhang et al. Reference Zhang, Xu, Liu and Dang2020). When the heating parameter is gradually increased, heat transfer first deteriorates, then recovers, and finally can approach as large values as for downward flow (Zhang et al. Reference Zhang, Xu, Liu and Dang2020). Interestingly, shear turbulence is gradually suppressed and even laminarised at lower Reynolds numbers (Bae et al. Reference Bae, Yoo and Choi2005; Chu, Laurien & McEligot Reference Chu, Laurien and McEligot2016; He, He & Seddighi Reference He, He and Seddighi2016; He et al. Reference He, Tian, Jiang and He2021), then the flow enters a convective turbulence state when the buoyancy exceeds a critical intensity (Su & Chung Reference Su and Chung2000; Marensi, He & Willis Reference Marensi, He and Willis2021).
Research on the phenomenon of laminarisation in mixed convection can be traced at least as far back as Hall & Jackson (Reference Hall and Jackson1969), who provided a theoretical explanation of this phenomenon, suggesting that reduced shear stress in the buffer layer leads to a reduction or even elimination of turbulence production. This interpretation received wide acceptance. More recently, He et al. (Reference He, He and Seddighi2016) modelled the buoyancy with a radially dependent body force added to isothermal flow, successfully reproducing the laminarisation phenomenon. They found that the body force causes little change to the key characteristics of turbulence, and proposed that laminarisation is caused by the reduction of the ‘apparent Reynolds number’, which is calculated based only on the pressure force of the flow (i.e. excluding the contribution of the body force). In flows at supercritical pressure, the laminarisation and deterioration of heat transfer are reproduced successfully by Bae et al. (Reference Bae, Yoo and Choi2005, Reference Bae, Yoo, Choi and McEligot2006). He et al. (Reference He, Tian, Jiang and He2021) researched the laminarisation phenomenon in flows at supercritical pressure, and established a unified explanation for the laminarisation mechanism. It is thought to be due to the variations of thermophysical properties, buoyancy and inertia, the latter of which plays a significant role in a developing flow.
Meanwhile, in isothermal flow, the laminarisation phenomenon has been observed when the base velocity profile is flattened (Hof et al. Reference Hof, De Lozar, Avila, Tu and Schneider2010; Kühnen et al. Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018; Marensi, Willis & Kerswell Reference Marensi, Willis and Kerswell2019). (The effect of buoyancy on the base velocity profile can be similar.) Kühnen et al. (Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018) proposed that the laminarisation is caused by reduced transient growth of small perturbations when the velocity profile is flattened. Marensi et al. (Reference Marensi, Willis and Kerswell2019) considered the nonlinear stability of finite perturbations for a series of flattened base velocity profiles in (isothermal) pipe flow, and found that flattening enhances the nonlinear stability of the laminar flow. Recently, machine learning has been used to explore laminarisation events in a reduced model of wall-bounded shear flows, and it has been suggested that the collapse of turbulence is connected with the suppression of streak instabilities (Lellep et al. Reference Lellep, Prexl, Eckhardt and Linkmann2022).
There have been many studies of mixed convection, both experimental (Celataa et al. Reference Celataa, Dannibale, Chiaradia and Cumo1998; Wang et al. Reference Wang, Li, Yu and Chen2011; Jackson Reference Jackson2013; Zhang et al. Reference Zhang, Xu, Liu and Dang2020) and numerical (You, Yoo & Choi Reference You, Yoo and Choi2003; Poskas, Poskas & Gediminskas Reference Poskas, Poskas and Gediminskas2012; Yoo Reference Yoo2013; Zhao et al. Reference Zhao, Zhu, Ge, Liu and Li2018), but most of these works focus on statistical properties, e.g. the Nusselt number, Reynolds stress, mean velocity profile and the mean temperature profile. Few works pay attention to the dynamical characteristics of the flow, such as the transition between flow types and the maintenance of convective turbulence. In particular, the transition mechanisms for mixed convection appear to be quite different from those of isothermal flow. Experimental results (Hanratty, Rosen & Kabel Reference Hanratty, Rosen and Kabel1958; Scheele, Rosen & Hanratty Reference Scheele, Rosen and Hanratty1960; Ackerman Reference Ackerman1970) have shown that a heated vertical pipe flow can go through a flow transition at Reynolds number as low as 
$30$. Scheele & Hanratty (Reference Scheele and Hanratty1962) found that the vertical heated upward flow in a pipe first becomes unstable when the velocity profiles develop an inflectional point. The flow develops regular and periodic motion after transition at a rather low Reynolds number. Similar patterns have been also observed by Kemeny & Somers (Reference Kemeny and Somers1962). Yao (Reference Yao1987) confirmed the experimental observations that the flow in a heated vertical pipe is supercritically unstable, and found that the bifurcated new equilibrium laminar flow is likely to be a double spiral flow. A weakly nonlinear instability analysis by Rogers & Yao (Reference Rogers and Yao1993) revealed that the heated upward flow is supercritically unstable while the heated downward flow is potentially subcritically unstable. Su & Chung (Reference Su and Chung2000) systematically researched the linear stability of mixed convection for upward and downward flow in a vertical pipe with constant heat flux, to explain multiple flow states that appeared in the experimental and numerical results. The calculation found that the first azimuthal mode is always the most unstable. The Rayleigh–Taylor instability is operative when the Reynolds number is extremely low, while the opposed-thermal-buoyant instability is dominant when the Reynolds number is higher. The transition of a low Reynolds number mixed convection flow in a vertical channel has been investigated by Chen & Chung (Reference Chen and Chung2002) for K-type (Klebanoff, Tidstrom & Sargent Reference Klebanoff, Tidstrom and Sargent1962) and H-type (Herbert Reference Herbert1983) disturbances. It was found that the flow field bifurcates to a new quasi-steady nonlinear state after the initial transient period, for both types of perturbation. Khandelwal & Bera (Reference Khandelwal and Bera2015) developed a weakly nonlinear stability theory in terms of a Landau equation to analyse the nonlinear saturation of stably stratified non-isothermal Poiseuille flow in a vertical channel with respect to different fluids, i.e. mercury, gases, liquids and heavy oils. A substantial enhancement in heat transfer rate was found for liquids and heavy oils from the basic state beyond the critical Rayleigh number. Recently, Marensi et al. (Reference Marensi, He and Willis2021) also performed linear stability analysis for a vertical heated pipe flow, using the parameter 
$C$ to measure the buoyancy force relative to the force that drives the laminar isothermal shear flow. Their results show that the flow is not always linearly unstable at a strong buoyancy condition, and revealed that the heat transfer deterioration is caused by the suppression of rolls in the flow. This finding is in line with the study of Lellep et al. (Reference Lellep, Prexl, Eckhardt and Linkmann2022), as the suppression of rolls would occur when streak instabilities are suppressed.
Some works have considered the linear and weakly nonlinear stability of a vertical heated flow, and have shown the instability of convective turbulence. However, the transition involving shear-driven turbulence is fundamentally nonlinear, and multiple flow states can arise at the same parameter values. Further understanding is of great significance for heat transfer prediction and optimisation in engineering applications, and there remain open questions to be addressed. How does the shear turbulence gradually disappear? Is laminarisation similar to that in isothermal flow? How is the convective turbulence triggered and sustained?
To tackle these questions, we employ nonlinear non-modal stability analysis (Pringle & Kerswell Reference Pringle and Kerswell2010; Cherubini et al. Reference Cherubini, De Palma, Robinet and Bottaro2011; Pringle, Willis & Kerswell Reference Pringle, Willis and Kerswell2012; Cherubini & De Palma Reference Cherubini and De Palma2013; Marensi et al. Reference Marensi, Willis and Kerswell2019) to seek the most efficient perturbation to trigger transition from the laminar state in a vertical heated pipe. The optimised perturbation is called the nonlinear optimal (NLOP), and the critical (i.e. lowest energy) NLOP that triggers turbulence is called the ‘minimal seed’ in isothermal flow. The magnitude and structure of the minimal seed reflects the nonlinear stability properties of the flow (Marensi et al. Reference Marensi, Willis and Kerswell2019). The Lagrange multiplier method has been used widely in fluid mechanics, including the heat transfer optimisation approach proposed by Guo, Zhu & Liang (Reference Guo, Zhu and Liang2007) and Motoki, Kawahara & Shimizu (Reference Motoki, Kawahara and Shimizu2018). We do not optimise heat transfer directly here, however, but seek the optimal flow perturbation that leads to transition away from the low heat-transfer state of laminar flow.
The plan of the paper is as follows. In § 2, we present the new model for direct numerical simulations of vertical heated pipe flow, and methods for the linear stability and nonlinear non-modal stability analysis. In § 3, we first show the results of direct numerical simulations, then present the results of linear stability and nonlinear non-modal stability analysis. Next, we reveal the separate transition processes starting from the NLOP and unstable eigenfunctions, as well as the self-sustaining process in convective turbulence. Finally, the paper concludes with a summary in § 4.
2. Formulation
2.1. Simulation of heated pipe flow
The upward flow through a vertical heated pipe flow is considered in this work. Like several models, we assume that there is a background temperature gradient along the axis of the pipe – e.g. Yao (Reference Yao1987), Chen & Chung (Reference Chen and Chung1996, Reference Chen and Chung2002) and Khandelwal & Bera (Reference Khandelwal and Bera2015). Different from many models, however, we suppose that the temperature gradient may vary in time due to changes in the way the fluid flows. We fix the difference between the bulk temperature of the fluid and that of the wall, aiming to see directly how the buoyancy parameter affects the flow pattern, and hence the heat flux at the wall and the temperature gradient. Our code is a modification to that of Marensi et al. (Reference Marensi, He and Willis2021), where the function of the temperature gradient was replaced by a spatially uniform heat sink. That assumption had the advantage of providing a simple analytic expression for the laminar state; here, it is computed numerically.
Let 
$\boldsymbol {x}=(r,\phi,z)$ denote cylindrical coordinates. We decompose the total temperature as
\begin{equation} T_{tot}(\boldsymbol{x},t)=T_w(z,t) + T(\boldsymbol{x},t) - 2 T_b, \end{equation}
with wall temperature 
$T_w=a_{tot}(t)\,z+b$, where 
$a_{tot}(t)$ is the time-dependent axial temperature gradient, and 
$b$ is a constant reference temperature. Axial periodicity over a distance 
$L$ is assumed for the temperature fluctuation field 
$T(\boldsymbol {x},t)$ and velocity field 
$\boldsymbol {u}_{tot}$. Let 
$T_b=\langle T\rangle$, where the angle brackets denote the spatial average. The factor 
$-2T_b$ has been inserted in (2.1) so that the temperature fluctuations due to the heating from the wall are positive and largest at the wall: evaluating (2.1) at the pipe radius 
$r=R$ gives 
$T|_{r=R}=2T_b$.
We take 
$2T_b$ as our temperature scale for non-dimensionalisation, and assume that local fluctuations in 
$T$ are more rapid than changes in the gradient 
$a_{tot}(t)$. For the length and velocity scales, we non-dimensionalise with the radius 
$R$ and twice the bulk flow speed 
$2U_b$, the latter of which coincides with the isothermal laminar centreline speed. Throughout the rest of the text, dimensionless variables and equations are presented, except in the definition of dimensionless parameters.
Temperature fluctuations 
$T$ satisfy
\begin{equation} \frac{\partial T }{\partial t}+({\boldsymbol{u}_{tot}}\boldsymbol{\cdot}\boldsymbol{\nabla}) T =\frac{1}{Re\,Pr}\,{\nabla} ^{2}T -{\boldsymbol{u}_{tot}}\boldsymbol{\cdot}{\hat{\boldsymbol{z}}}\,a_{tot}(t), \end{equation}
with boundary condition 
$T=1$ at the wall. A fixed bulk temperature 
$\langle T\rangle =1/2$ is maintained through adjustments in the gradient 
$a_{tot}(t)$, which is determined at each instant via the spatial average of the equation (see (2.11a,b) below). The dimensionless parameters are the Reynolds and Prandtl numbers 
$Re=2U_bR/\nu$ and 
$Pr=\nu /\kappa$, where 
$\nu$ and 
$\kappa$ are the kinematic viscosity and thermal diffusivity, respectively. Under the Boussinesq approximation (Turner & Turner Reference Turner and Turner1979), the Navier–Stokes equations are
\begin{equation} \frac{\partial\boldsymbol{u}_{tot}}{\partial t}+(\boldsymbol{u}_{tot}\boldsymbol{\cdot} \boldsymbol{\nabla}) {\boldsymbol{u}_{tot}}={-}\boldsymbol{\nabla} p_{tot} + \frac{1}{Re}\,{\nabla}^{2} {\boldsymbol{u}_{tot}}+\frac{4}{Re}\,(1+\beta'(t)+C T )\hat{\boldsymbol{z}} , \end{equation}with continuity equation
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_{tot}=0 , \end{equation}
and no-slip condition 
$\boldsymbol {u}_{tot}=\boldsymbol {0}$ at the wall. Here, 
$\beta '(t)$ is the excess pressure fraction, relative to isothermal laminar flow, required to maintain the fixed dimensionless mass flux, 
$\langle {\boldsymbol {u}_{tot}}\boldsymbol {\cdot }{\hat {\boldsymbol {z}}} \rangle =1/2$. The dimensionless parameter
\begin{equation} C=\frac{Gr}{16\,Re} \end{equation}
measures the buoyancy force relative to the force that drives laminar isothermal shear flow, where 
${Gr}={\gamma g(T|_{r=R}-T_b)(2R)^{3}}/{\nu ^{2}}$ is the Grashof number, 
$\gamma$ is the coefficient of volume expansion, and 
$g$ is gravitational acceleration. For further details, see Marensi et al. (Reference Marensi, He and Willis2021).
As the focus of this study is on the dynamics of perturbations from the laminar solution, we decompose the variables, 
${\boldsymbol {u}_{tot}}=\boldsymbol {u}_0+\boldsymbol {u}$, 
$p_{tot}=p_0+p$, 
$1+\beta '(t)=1+\beta _{0}+\beta (t)$, 
$T=\varTheta _0+\varTheta$, 
$a_{tot}(t)=a_{0}+a(t)$, where subscript 
$0$ tags the laminar solution, 
$\boldsymbol {u}=(u_r,u_{\phi },u_z)$. The laminar velocity 
$\boldsymbol {u}_0=u_0(r)\,\hat {\boldsymbol {z}}$ satisfies
\begin{equation} -\boldsymbol{\nabla} p_{0} + \frac{1}{Re}\,{\nabla}^{2} {\boldsymbol{u}_{0}}+\frac{4}{Re}\,(1+\beta_0+C \varTheta_0)\hat{\boldsymbol{z}}=0 . \end{equation}Subtracting (2.6) from (2.3), we obtain
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t}+u_0\,\frac{\partial \boldsymbol{u}}{\partial z}+u_{r}\,\frac{{\rm d}u_0}{{\rm d}r}\, \hat{\boldsymbol{z}}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} ={-}\boldsymbol{\nabla} p + \frac{1}{Re}\,{\nabla}^{2}\boldsymbol{u}+\frac{4}{Re}\,(C \varTheta+\beta(t))\hat{\boldsymbol{z}}. \end{equation}
Taking the spatial average of the 
$z$-components of (2.6) and (2.7) gives equations for 
$\beta _0$ and 
$\beta (t)$ that fix 
$\langle u_z \rangle = 0$:
\begin{equation} \beta_0 ={-}\frac{1}{2} \left( \left.\frac{\partial u_0}{\partial r}\right|_{r=1}+C\right), \quad \beta(t) ={-}\frac{1}{2} \left.\frac{\partial (u_z)_{00}}{\partial r}\right|_{r=1}, \end{equation}
where 
$(\cdot )_{00}$ corresponds to averaging over 
$\phi$ and 
$z$. The laminar temperature profile 
$\varTheta _0(r)$ satisfies
\begin{equation} \frac{1}{Re\,Pr}\,{\nabla}^{2}\varTheta_0 =\boldsymbol{u}_{0}\boldsymbol{\cdot}{\hat{\boldsymbol{z}}}a_{0} . \end{equation}Subtracting (2.9) from (2.2), we obtain
\begin{equation} \frac{\partial \varTheta }{\partial t}+u_0\,\frac{\partial \varTheta }{\partial z}+u_r\,\frac{{\rm d}\varTheta_{0}} {{\rm d}r} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \varTheta =\frac{1}{Re\,Pr}\,{\nabla}^{2}\varTheta -u_z a_{0}-(u_0+u_z)\,a(t). \end{equation}
Taking the spatial averages of (2.9) and (2.10) gives equations that determine 
$a_0$ and the value of 
$a(t)$ required to fix 
$\langle \varTheta \rangle =0$:
\begin{equation} a_0 = \frac{4}{Re\,Pr} \left.\partial_{r}\varTheta_{0}\right|_{r=1},\quad a(t)= \frac{4}{Re\,Pr} \left.\partial_{r}(\varTheta)_{00}\right|_{r=1} . \end{equation}2.2. Linear stability analysis
Arnoldi iteration is employed to calculate the leading eigenvalues and eigenfunctions of the laminar solution using our simulation code. Considering a small perturbation 
$\boldsymbol {u}$ to the laminar solution 
$\boldsymbol {u}_0$, so that the nonlinear terms of (2.7) may be ignored, the linearised system may be written in the form
\begin{equation} \partial _{t}\boldsymbol{u}=A(u_0)\,\boldsymbol{u}. \end{equation}
Integrating over a period 
$\mathcal {T}$, eigenfunctions of (2.12) with growth rate 
$\sigma$ satisfy the exponentiated eigenvalue problem
\begin{equation} \boldsymbol{u}(\mathcal{T}) = \boldsymbol{u}(0) + \int_0^{\mathcal{T}} A(u_0)\, \boldsymbol{u}(t) \, \mathrm{d}t = B(u_0)\,\boldsymbol{u}(0) = \mathrm{e}^{\sigma \mathcal{T}}\boldsymbol{u}(0) . \end{equation}
The Arnoldi method requires only the result of the calculations of multiplying 
$B$ with given 
$\boldsymbol {u}$, i.e. the result of time integration of 
$\boldsymbol {u}$ over the period 
$\mathcal {T}$. Given a starting 
$\boldsymbol {u}$, the method seeks eigenvectors in the Krylov subspace 
$\mathcal {K}=\mathrm {span} \{\boldsymbol {u}, B\boldsymbol {u}, B^2\boldsymbol {u},\ldots \}$, using Gram–Schmidt orthogonalisation to improve the numerical suitability of this basis set. Using 
$f(\boldsymbol {u})$ to denote the result of time integration of a perturbation 
$\boldsymbol {u}$ from the laminar state, we may make the approximation
\begin{equation} B(u_0)\, \boldsymbol{u}\approx \frac{1}{\epsilon } \left\{ \,f(\boldsymbol{0}+\epsilon \boldsymbol{u})-f(\boldsymbol{0}) \right\} =\frac{1}{\epsilon}\,f(\epsilon\boldsymbol{u}) , \end{equation}
for some small value 
$\epsilon$. We take 
$\epsilon$ such that 
$\epsilon \,\|\boldsymbol {u}\|=10^{-6}\,\|\boldsymbol {u}_0\|$, i.e. the perturbation is 6 orders of magnitude smaller than the laminar state, and nonlinear terms are 12 orders smaller. The choice of 
$\mathcal {T}$ does not affect evaluation of the real part of 
$\sigma$, which determines stability, but a large 
$\mathcal {T}$ can cause aliasing of its complex part. The period 
$\mathcal {T}$ is therefore chosen to be small, but not too small, as the number of Arnoldi iterations required is usually inversely proportional to 
$\mathcal {T}$, and many iterations would require storage of a large 
$\mathcal {K}$. With the aim of keeping the calculation to approximately 50 Arnoldi iterations, typically 
$\mathcal {T}\approx 10$.
2.3. Nonlinear non-modal stability analysis
When the linear system is stable, we must appeal to nonlinear dynamics to interpret transition. The NLOP (Pringle & Kerswell Reference Pringle and Kerswell2010; Pringle et al. Reference Pringle, Willis and Kerswell2012) is the optimal perturbation that achieves the largest energy growth, and when it is just large enough to trigger transition, it is called the minimal seed.
The NLOP can be found by the Lagrange multiplier technique; for details, see Pringle et al. (Reference Pringle, Willis and Kerswell2012). To reduce complexity, we consider only the optimal velocity perturbation, and set the initial temperature perturbation to zero. In this case, the Lagrangian is defined by
\begin{align} \mathcal{L} &= \big\langle\tfrac{1}{2}\,(\boldsymbol{u}(x,\mathcal{T}))^{2}\big\rangle - \lambda_0 \big(\tfrac{1}{2}\,\langle(\boldsymbol{u}(x,0))^{2}\rangle-E_{0}\big)- \int_{0}^{\mathcal{T}}\langle \boldsymbol{v} \boldsymbol{\cdot} \mathrm{NS}(\boldsymbol{u})\rangle \, \mathrm{d}t\nonumber\\ &\quad -\int_{0}^{\mathcal{T}}\langle\varPi \, \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}\rangle \, \mathrm{d}t -\int_{0}^{\mathcal{T}}\langle{\rm \pi} \, \mathrm{Tem}(\varTheta)\rangle \, \mathrm{d}t\nonumber\\ &\quad - \int_{0}^{\mathcal{T}} \varGamma \langle \boldsymbol{u}\boldsymbol{\cdot} \hat{\boldsymbol{z}} \rangle\, \mathrm{d}t-\int_{0}^{\mathcal{T}} Q \langle \varTheta \rangle\, \mathrm{d}t , \end{align}
where 
$\lambda _0$, 
$\varPi (\boldsymbol {x},t)$, 
${\rm \pi} (\boldsymbol {x},t)$, 
$\varGamma (t)$, 
$Q(t)$ and 
$\boldsymbol {v}(\boldsymbol {x},t)=(v_r,v_{\phi },v_z)$ are Lagrange multipliers. The first term, the perturbation energy at final time 
$t=\mathcal {T}$, is the objective function to be maximised. The second term is the constraint of fixed amplitude for the initial perturbation. Then the velocity perturbation 
$\boldsymbol {u}$ is constrained to satisfy the Navier–Stokes equation 
$\mathrm {NS}(\boldsymbol {u})$ and the continuity equation from 
$t=0$ to 
$t=\mathcal {T}$, while the temperature perturbation satisfies the temperature equation 
$\mathrm {Tem}(\varTheta )$ from 
$t=0$ to 
$t=\mathcal {T}$. The velocity and temperature must also satisfy fixed mass flux and fixed heat flux, represented by the last two terms. Taking variations of 
$\mathcal {L}$ with respect to each variable, and setting them equal to zero, we obtain the following set of Euler–Lagrange equations. The adjoint Navier–Stokes, temperature equation and continuity equations are
\begin{gather} \frac{\delta \mathcal{L}} {\delta \boldsymbol{u}}=\frac{\partial \boldsymbol{v}}{\partial t}+u_0\,\frac{\partial \boldsymbol{v}}{\partial z}-v_{z}{u'_0}{\hat{\boldsymbol{r}}}+ \boldsymbol{\nabla}\times(\boldsymbol{v}\times \boldsymbol{u})- \boldsymbol{v}\times\boldsymbol{\nabla} \times \boldsymbol{u} +\boldsymbol{\nabla} \varPi\nonumber\\ +\frac{1}{Re}\,{\nabla}^{2}\boldsymbol{v}-{\rm \pi}\,\varTheta'_{0}\hat{\boldsymbol{r}}-{\rm \pi}\boldsymbol{\nabla}\varTheta-{\rm \pi} (a(t)+a_0)\hat{\boldsymbol{z}}-\varGamma \hat{\boldsymbol{z}}=0 , \end{gather}
\begin{gather} \frac{\delta \mathcal{L}}{\delta \varTheta}=\frac{\partial {\rm \pi}}{\partial t}+u_0\,\frac{\partial {\rm \pi}}{\partial z}+\frac{4}{Re}\,v_z C+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}{\rm \pi} + \frac{1}{Re\, Pr}\,{\nabla}^2{\rm \pi}-Q=0 , \end{gather}
\begin{gather} \frac{\delta \mathcal{L}}{\delta p}=\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v} =0 , \end{gather}where primes indicate the radial derivative. The compatibility conditions are given by
$$\begin{gather} \frac{\delta \mathcal{L}}{\delta\boldsymbol{u}(\boldsymbol{x},\mathcal{T})}=\boldsymbol{u}(\boldsymbol{x},\mathcal{T})-\boldsymbol{v}(\boldsymbol{x},\mathcal{T})=0, \end{gather}$$
$$\begin{gather}\frac{\delta \mathcal{L}}{\delta \varTheta(\boldsymbol{x},\mathcal{T})}={-}{\rm \pi}(\boldsymbol{x},\mathcal{T})=0, \end{gather}$$and the optimality condition is
\begin{equation} \frac{\delta \mathcal{L}}{\delta {\boldsymbol{u}(\boldsymbol{x},0)}}={-}\lambda_0\, \boldsymbol{u}(\boldsymbol{x},0)+\boldsymbol{v}(\boldsymbol{x},0)=0. \end{equation}
Optimisation is performed iteratively from a starting velocity field 
$\boldsymbol {u}(\boldsymbol {x},0)^{(0)}$, where the superscript indicates the iteration. Given a field 
$\boldsymbol {u}(\boldsymbol {x},0)^{(j)}$, the Navier–Stokes equations are integrated to time 
$\mathcal {T}$. The compatibility conditions (2.19) and (2.20) provide end conditions for the backwards integration of (2.16) and (2.17) to time 
$0$. Like 
$p$ and 
$\beta$, the Lagrange multipliers 
$\varPi$ and 
$\varGamma$ are used to ensure that 
$\boldsymbol {\nu }(\boldsymbol {x},t)$ is also divergence-free and has zero flux. Here, 
$Q$ is used to fix 
$\langle {\rm \pi}\rangle =0$. Finally, once 
$\boldsymbol {\nu }(\boldsymbol {x},0)$ has been calculated, the optimality condition (2.21) may be evaluated. As the optimality condition is not satisfied automatically, 
$\boldsymbol {u}(\boldsymbol {x},0)$ is moved in the direction of 
$\delta \mathcal {L}/\delta \boldsymbol {u}(\boldsymbol {x},0)$ towards achieving a maximum for 
$\mathcal {L}$, where 
$\delta \mathcal {L}/\delta \boldsymbol {u}(\boldsymbol {x},0)$ should vanish. The update to 
$\boldsymbol {u}(\boldsymbol {x},0)^{(j)}$ each iteration is
\begin{equation} \boldsymbol{u}(\boldsymbol{x},0)^{(j+1)}=\boldsymbol{u}(\boldsymbol{x},0)^{(j)}+{\epsilon_0}\,\frac{\delta \mathcal{L}}{\delta \boldsymbol{u}(\boldsymbol{x},0)^{(j)}} , \end{equation}
where 
$\epsilon _0$ is a small value, controlled using a procedure described in Pringle et al. (Reference Pringle, Willis and Kerswell2012), and 
$\lambda _0$ is adjusted to set 
$\langle [\boldsymbol {u}(\boldsymbol {x},0)^{(j+1)}]^2\rangle =2\,E_0$.
2.4. Time integration code
The calculations are carried out by the open-source code openpipeflow.org (Willis Reference Willis2017). Variables are discretised on the domain 
$\{ r,\phi,z \}=[0,1]\times [0,2{\rm \pi} ]\times [0,2{\rm \pi} /k_0]$, where 
$k_0=2{\rm \pi} /L$, using Fourier decomposition in the azimuthal and streamwise directions and finite difference in the radial direction,
\begin{equation} \left\{\boldsymbol{u},p,\varTheta\right\}(r_s,\phi,z)=\sum_{k<|K|} \sum_{m<|M|}\left\{\boldsymbol{u},p,\varTheta \right\}_{skm}\exp({{\rm i}(k_0kz+m\phi)}) , \end{equation}
where 
$s=1,\ldots,S$, and the radial points are clustered near the wall. Temporal discretisation is via a second-order predictor–corrector scheme, with an Euler predictor and a Crank–Nicolson corrector applied to the nonlinear terms. The laminar solution is quickly calculated by eliminating azimuthal and axial variations using a resolution 
$S=64$, 
$M=1$, 
$K=1$. For nonlinear calculations, as the adjoint optimisation requires forward and backward time integrations for each iteration, it is computationally expensive. To keep calculations manageable, a Reynolds number 
$Re=3000$ is adopted for the nonlinear non-modal stability analysis with a domain of length 
$L=5D$. We use 
$S=64$, 
$M=48$, 
$K=42$, and a time step 
$\Delta t=0.01$.
3. Results
3.1. Flow regimes
Figure 1 shows how the laminar velocity and temperature profile changes with 
$C$. It is observed that the velocity profile becomes flattened, then is well known to become ‘M-shaped’ at larger values of 
$C$. The temperature profile also changes, but not so strongly. The boundary temperature gradient of the laminar state is shown in figure 2(a). The laminar temperature gradient increases slightly with 
$C$, and decreases substantially with 
$Re$, as the higher flow rate carries heat out of the system more rapidly. The Nusselt number is defined as
\begin{equation} {Nu}=\frac{2R\, q_w}{\lambda ({T|_{r=R}-T_b})}, \end{equation}
where 
$\lambda$ is the thermal conductivity of the fluid, and 
$q_w$ is net wall heat flux. The laminar Nusselt number increases with a larger value of 
$C$ (see figure 2b), but it is independent of the Reynolds number (Su & Chung Reference Su and Chung2000).
Figure 1. (a) The laminar velocity profile. (b) The laminar temperature profile. Laminar profiles are independent of Reynolds number.
Figure 2. (a) Temperature gradient of laminar flow 
$a_0$ as a function of 
$C$ for several 
$Re$. (b) Nusselt number as a function of 
$C$.
As discussed in the Introduction, there are three typical flow states: laminar flow, shear turbulence and convective turbulence (convection-driven turbulence). We have performed a suite of simulations at different Reynolds numbers 
$Re$ and buoyancy parameters 
$C$, and the observed flow regime is shown in figure 3(a). The diagram is consistent with figure 7 of Marensi et al. (Reference Marensi, He and Willis2021). Close to the boundary, where the laminar or convective turbulence regimes meet shear turbulence, the flow can be found in either state. Nevertheless, the three typical flow regimes can be identified clearly. In the isothermal case (
$C=0$), shear turbulence appears for 
$Re\gtrsim 2000$ (Avila et al. Reference Avila, Moxey, De Lozar, Avila, Barkley and Hof2011; Avila, Barkley & Hof Reference Avila, Barkley and Hof2023), but increasing the effect of buoyancy, at lower Reynolds numbers, the flow first laminarises then transitions to convective turbulence. At a higher Reynolds number, the laminarisation regime becomes narrower and finally disappears. Convective turbulence first appears for 
$C\gtrsim 4$ at lower Reynolds number, as reported by Scheele et al. (Reference Scheele, Rosen and Hanratty1960), Scheele & Hanratty (Reference Scheele and Hanratty1962) and Yao (Reference Yao1987). But at high Reynolds number, the shear-driven state persists to larger 
$C$ and crosses the critical 
$C$ observed at low Reynolds number, setting up a direct transition from the shear-driven turbulence to convective turbulence as 
$C$ is increased. The flow states at 
$Re$ larger than 5300 have not been computed in the present work, but this direct transition is known to occur at much larger 
$Re$, and a collapse in experimental data is observed when presented in terms of the buoyancy number 
$Bo= 8\times 10^4\, (8\,Nu\,Gr)/(Re^{3.425}\,Pr^{0.8})$, where 
$Gr = 16\,Re\,C$, as seen in figure 3(b).
Figure 3. (a) Regions of laminar flow (L), shear turbulence (S) and convective turbulence (C); SL and SC indicate that the flow may be in either of the two states. (b) Change in heat flux, normalised by that for the isothermal state (
$C\to 0$), as a function of 
$Bo = 8\times 10^4\,(8\,Nu\,Gr)/(Re^{3.425}\,Pr^{0.8})$. Present data from simulations at 
$Re = 5300$, 
$Pr = 0.7$ and various 
$Gr = 16\,Re\,C$. The upper and lower branches correspond to shear turbulence and convective turbulence, respectively.
Numerical results are validated quantitatively by comparison with experiments in figure 3(b), which also includes the numerical results of Marensi et al. (Reference Marensi, He and Willis2021). The strong heat flux for smaller values of buoyancy number 
$Bo$ is brought about by shear turbulence, but deteriorates as 
$Bo$ is increased when buoyancy suppresses the shear turbulence. When 
$Bo$ exceeds a certain value, the flow switches from shear turbulence to convective turbulence, accompanied by a sudden drop in heat transfer. However, the heat transfer eventually recovers with further increase of 
$Bo$. The captured heat transfer features are consistent with the results reported by Bae et al. (Reference Bae, Yoo and Choi2005) and Yoo (Reference Yoo2013). In the present model, the velocity field is allowed to affect the heat flux at the wall; the model and code are a minor update to that of Marensi et al. (Reference Marensi, He and Willis2021), which assumed a constant boundary temperature with a spatially uniform heat sink term. The uniform sink has the advantage that there is a simple analytic expression for the laminar base flow, but the axial temperature gradient that leads to the spatially dependent heat sink term is expected to be closer to the real system. Compared with Marensi et al. (Reference Marensi, He and Willis2021), it is found that the present model provides a small improvement in agreement with experimental data in the shear turbulence regime, while in the convective turbulence regime, the two models give similar results.
Starting from isothermal flow (
$C=0$) at 
$Re=3000$, figure 4 shows the effect of switching on the buoyancy for 
$C=1,3,5$. When the buoyancy is weak, at 
$C=1$, in figure 4(a), we find that all three components maintain a high amplitude, and the flow is essentially unchanged from the isothermal case. For 
$C=3$, the velocity perturbation energy decays quickly, indicating the occurrence of laminarisation. At 
$C=5$, initially a similar energy drop is also observed, but it stops at a much lower energy level, and fluctuates with a low frequency, which is typical of the convective turbulence state. Heat transfer for the laminar and convective turbulence state is severely reduced, seen in figure 4(b), where the Nusselt number drops to approximately half of that for shear turbulence. The convective turbulence state at these parameters has a Nusselt number that is not much greater than for the laminar state.
Figure 4. Simulation at 
$Re=3000$ starting from the same initial condition at 
$C=1, 3, 5$. (a) Velocity perturbation energy. The solid, dashed and dotted lines represent the energy of three components of velocity, denoted as 
$E(u_r)$, 
$E(u_{\phi })$ and 
$E(u_z)$, respectively. (b) Nusselt number.
Isosurfaces of the streamwise vorticity are plotted in figure 5 for shear turbulence (
$C=1$) and convective turbulence (
$C=5$). Shear turbulence has complicated streamwise vortices that fill the whole pipe, while the convective turbulence exhibits more organised vortex structures that are concentrated mainly around the middle of the pipe. The location of the structures can be seen more clearly in the cross-sections of figure 6. In shear turbulence, there are abundant near-wall streamwise vortices, which lift up low-speed streaks, and enhance the mixing of fluid in the pipe. In convective turbulence, the perturbations are found around the centre of the pipe, and the near-wall flow is almost steady. There are almost no near-wall rolls and streaks observed. Further visualisations of flow at larger 
$C$ show similar characteristics. The self-sustaining process of shear turbulence, reliant on near-wall vortices for the creation of streaks, has been destroyed in convective turbulence. The lack of streamwise vortices significantly weakens the interaction between fluid near the wall and centre, and heat transfer is therefore seriously damaged.
Figure 5. The isosurface of streamwise vorticity at (a) shear turbulence (
$C=1$) and (b) convective turbulence (
$C=5$) at 
$Re=3000$; red/yellow are 
$30\,\%$ of the min/max streamwise vorticity.
Figure 6. Cross-sections of (a) shear turbulence (
$C=1$) and (b) convective turbulence (
$C=5$) at 
$Re=3000$. The contours are coloured by streamwise total velocity, and the arrows represent cross-stream components. The largest arrow has magnitude 
$5.6\times 10^{-3}$ in (a), and 
$1.93\times 10^{-4}$ in (b).
3.2. Linear stability results
This subsection examines the change of linear stability of the laminar solution at different 
$C$ values, revealing how the buoyancy affects the dynamics near the laminar solution. Time stepping is combined with Arnoldi iteration to calculate the leading eigenvalues and eigenfunctions for small perturbations about the laminar state. A resolution 
$S=64$, 
$M=10$, 
$K=10$ is used. A finer radial resolution does not produce a noticeable difference in the results, and the low 
$K$ and 
$M$ are sufficient to pick the most unstable axial and azimuthal modes for the given domain (
$L=5D$). It is found that the most unstable mode is always of azimuthal wavenumber 
$m=1$ and is usually of wavenumber 
$k=1$ or 
$2$. These observations are consistent with results of Scheele et al. (Reference Scheele, Rosen and Hanratty1960), Su & Chung (Reference Su and Chung2000) and Marensi et al. (Reference Marensi, He and Willis2021). The real part of the eigenvalue of the unstable modes is shown in figure 7. The instability first appears near 
$C=4$, consistent with the appearance of convective turbulence, verifying that convective turbulence is caused by linear instability (Yao Reference Yao1987; Rogers & Yao Reference Rogers and Yao1993; Marensi et al. Reference Marensi, He and Willis2021). Two branches of unstable mode lead to the instability of the flow, apart from at the lowest 
$Re$. For both 
$k=1$ and 
$k=2$ cases, the first branch enters just before 
$C=4$ but stabilises by 
$C\approx 6$. A second branch then takes over for 
$C\gtrsim 6$. Comparing eigenvalues between the 
$k=1$ and 
$k=2$ cases, the first branch with 
$k=1$ leads at lower 
$C$, and the second branch with 
$k=2$ leads at larger 
$C$. Interestingly, the flow is not always linearly unstable for 
$C>4$, as there is a window of stability for 
$C\approx 6$ for a range 
$Re\approx 2000\unicode{x2013}4000$. This was also reported for the model of Marensi et al. (Reference Marensi, He and Willis2021). In particular, at 
$Re=3000$, the flow is linearly stable for 
$C=6$ – we will take advantage of this for the nonlinear stability analysis of § 3.3. Eigenfunctions of the two branches of unstable modes with wavenumber 
$k=2$, 
$m=1$ are shown in figure 8 with isosurfaces of streamwise velocity and streamwise vorticity. The eigenfunctions are consistent with the results of Khandelwal & Bera (Reference Khandelwal and Bera2015) and Marensi et al. (Reference Marensi, He and Willis2021). The perturbation is active mainly in the centre of the pipe, in good agreement with the observations of Su & Chung (Reference Su and Chung2000) and Marensi et al. (Reference Marensi, He and Willis2021).
Figure 7. Growth rate of perturbations of different streamwise wavenumbers, (a) 
$m=1$, 
$k=1$, and (b) 
$m=1$, 
$k=2$, versus buoyancy parameter 
$C$ at several Reynolds numbers.
Figure 8. Isosurfaces of eigenfunctions for wavenumber 
$k=2$, 
$m=1$ at 
$Re=3000$: (a,b) first branch at 
$C=4$; (c,d) second branch at 
$C=10$, for (a,c) streamwise velocity perturbations, (b,d) streamwise vorticity. (Red/yellow are 30 % of the min/max values.)
Linear stability analysis results have confirmed that convective turbulence is triggered by linear instability, thus it appears mainly at 
$C>4$, as shown by the magenta shadowed region in figure 3(a). However, at a high Reynolds number, the critical 
$C$ for transition is pushed back to larger 
$C$. This indicates competition between shear turbulence and convective turbulence. For the limit of the shear turbulence region shown in figure 3(a), various possible explanations have been proposed. He et al. (Reference He, He and Seddighi2016) linked the buoyancy-modified flow to a partner isothermal flow at a different (lower) ‘apparent Reynolds number’, based on an apparent friction velocity associated with only the pressure force of the flow (i.e. excluding the contribution of the body force). They found that the buoyancy can reduce the apparent Reynolds number, thereby suppressing turbulence. Kühnen et al. (Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018) suggested that it is the flattened base velocity profile (see figure 1) that reduces the transient growth of streaky perturbations, and thereby causes the laminarisation. Marensi et al. (Reference Marensi, Willis and Kerswell2019) did a nonlinear non-normal stability analysis of the flattened base velocity profile, suggesting that the nonlinear stability of the flattened velocity profile is enhanced, i.e. a larger amplitude perturbation is required to trigger turbulence. Lellep et al. (Reference Lellep, Prexl, Eckhardt and Linkmann2022) used a machine learning method and linked the collapse of turbulence to the suppression of streak instability. We hold a view similar to that of Lellep et al. (Reference Lellep, Prexl, Eckhardt and Linkmann2022) – through linear and nonlinear optimisation (Chu, Willis & Marensi’ unpublished observations), we have found that a body force that flattens the base velocity profile always produces streaks with a more stable shape. Further numerical experiments have been carried out, imposing a body force that flattens the base velocity profile at particular times during the relative periodic orbit, i.e. only during the time interval of the formation of streaks, or the interval for the regeneration of rolls. Only when targeting the formation of streaks was the self-sustaining process disrupted, supporting that the flattened profile suppresses turbulence by the stabilisation of streaks. With the increase of Reynolds number, the streaks formed by the flattened base velocity profile will recover the instability to sustain turbulence. Therefore, the laminarisation regime gradually disappears at a higher Reynolds number; see the green shadowed region in figure 3(a).
3.3. Nonlinear non-modal stability analysis
Isothermal pipe flow is linearly stable, and nonlinear dynamics must be considered when examining transition. As 
$C$ is increased, we have just seen that linear instability arises, but also that there is a range of linear stability even in the convective regime. In this subsection, we consider a perturbation of given finite amplitude 
$A_0$ that grows most over a time 
$\mathcal {T}$. This perturbation is the NLOP. When 
$A_0$ is just large enough to trigger turbulence, the NLOP is known as the ‘minimal seed’. We seek to see how the NLOP is affected by buoyancy, and focus on 
$Re=3000$ for 
$C$ in the range 0–6. Nonlinear flows at 
$C=0,1,2$ are in the shear turbulence regime, while flow at 
$C=3$ returns to the laminar regime. Flows at 
$C=4,5,6$ are in the convective turbulence regime. Note, however, that only the laminar flows at 
$C=4$ and 
$5$ are linearly unstable, while at 
$C=6$ the flow is linearly stable yet well within the convective regime.
To calculate the largest growing perturbation to the laminar state, we use the method of § 2.3. Following Pringle et al. (Reference Pringle, Willis and Kerswell2012), 
$\mathcal {T}$ needs to be sufficiently large to produce a rapid jump in the objective function with respect to a relatively small increase in 
$E_0$. We first consider a target time 
$\mathcal {T}=300$, which is approximately twice the time at which peak linear growth can be observed at this Reynolds number. This 
$\mathcal {T}$ is found to be sufficient to isolate effects of buoyancy on the NLOP. Typical residuals 
$\delta \mathcal {L}/\delta \boldsymbol {u}(x,0)$ and energy growth 
$G=\langle \boldsymbol {u}(\mathcal {T})^2\rangle /\langle \boldsymbol {u}(0)^2\rangle$ during the optimisation are shown in figure 9(a). A clear drop of residual leads to final optimal energy growth, which means that the calculation is well converged, and starting each calculation from different random initial velocity fields produces the same final optimal. Usually, the calculation is stopped when the change in 
$G$ is smaller than 
$10\times 10^{-5}$.
Figure 9. (a) The residual 
$\delta \mathcal {L}/\delta \boldsymbol {u}(x, 0)$ (left-hand axis) and energy growth 
$G=\langle u(\mathcal {T})^2\rangle /\langle u(0)^2\rangle$ (right-hand axis) versus iterations in the case at 
$C=1$, 
$E_0=2.7\times 10^{-5}$, 
$Re=3000$. (b) The optimal energy growth 
$G$ changes with 
$E_0$ for 
$C=0\unicode{x2013}5$, 
$\mathcal {T}=300$, 
$Re=3000$.
We first present the results at 
$C=0\unicode{x2013}5$, where similar NLOPs are found, then a new NLOP at 
$C=6$ is presented separately. The optimal energy growth as a function of 
$E_0$ for several 
$C$ is plotted in figure 9(b). For each 
$C$ at small 
$E_0$, the growth 
$G$ is constant, which implies that the linear optimal has been found. As 
$E_0$ is increased, the NLOP takes over and 
$G$ increases. At slightly larger 
$E_0$ still, a clear jump in the energy growth is observed when shear turbulence is triggered. Note, however, that at 
$C=5$, the NLOP is not observed to take over for the given 
$\mathcal {T}=300$.
Figure 10 shows the time series of the energy growth of an NLOP (solid line) and minimal seed (dashed line) for 
$C=0\unicode{x2013}5$. The minimal seed triggers turbulence, while the perturbation eventually decays when it starts from the NLOP at a slightly lower 
$E_0$. (The true value of 
$E_0$ for the minimal seed could be further refined, but without a noticeable change in the velocity field. We will call the perturbation at the upper 
$E_0$ the minimal seed.) At a larger 
$C$, the maximum energy growth is smaller and the peak occurs at an earlier time. The reduced maximum energy is caused by the more flattened base velocity profile, which suppresses the nonlinear instability (Marensi et al. Reference Marensi, Willis and Kerswell2019). The maximum occurring at an earlier time is attributed to the smaller length scale of the streamwise vortices structure in the minimal seed, shown in figure 11(d).
Figure 10. The energy growth of the NLOP versus time at (a) 
$C=0\unicode{x2013}2$ and (b) 
$C=3\unicode{x2013}5$. In (b) for 
$C=5$, the dashed line is for an NLOP of the same type as for 
$C<5$, while for the solid line, the NLOP is close to the linear eigenfunction.
Figure 11. Isosurfaces of streamwise velocity perturbation and cross-sections at 
$z=5$ for the NLOP: (a,c) 
$C=0$, 
$E_0=1.9\times 10^{-5}$, and (b,d) 
$C=2$, 
$E_0=3.36\times 10^{-5}$, at 
$Re=3000$. Red/yellow isosurfaces are for 
$30\,\%$ of the min/max streamwise velocity. Arrows show the cross-stream components; the largest arrow has magnitude 
$1.55\times 10^{-4}$ in (c), and 
$2.39\times 10^{-4}$ in (d).
At 
$C=5$, the algorithm converges to either the solid green line or the dashed green line in figure 10(b), depending on the initial guess, although the solid green case appears more frequently. Note that the two optimals have similar energy at the target time 
$\mathcal {T}$, leading to the possibility of local optimals of the Lagrangian. (Local optimals with similar 
$\mathcal {L}$ have also been observed by Motoki et al. (Reference Motoki, Kawahara and Shimizu2018) for maximal heat transfer in steady plane Couette flow.) The solid green line includes a period of exponential growth, indicating that the perturbation approaches the unstable eigenfunction of the unstable laminar state. The initial state for the dashed line is of the same type as the NLOP calculated for 
$C<5$.
The structures of the minimal seeds at 
$C=0,2$ are shown in figure 11. They are essentially similar, localised in both the azimuthal and streamwise directions, consistent with the isothermal NLOP (Pringle & Kerswell Reference Pringle and Kerswell2010; Pringle et al. Reference Pringle, Willis and Kerswell2012). But compared with isothermal pipe flow 
$C=0$, the minimal seed at 
$C=2$ is more localised and ‘thinner’, being located closer to the wall; see figures 11(c,d). A thinner minimal seed was also reported by Marensi et al. (Reference Marensi, Willis and Kerswell2019) for a (prescribed) flattened base velocity profile. Here, the flattened profile is due to the buoyancy. The effective radial interval for the lift-up mechanism is narrower, and consequently the optimal rolls gradually become thinner. Structures of the NLOP at 
$C=3,4$ (not shown) are similar to figure 11(
$d$) but thinner still. The optimal perturbation at 
$C=5$, which is close to the unstable eigenfunction, is shown in figure 12(b), accompanied by the same branch of optimal perturbation at 
$C=4$ in figure 12(a), found by using a longer target time 
$\mathcal {T}=400$ (to enable slower linear growth to compete with the NLOP). This type of optimal perturbation is almost unchanged at small 
$E_0$, suggesting linearity. Its structures are initially also close to the wall, but after an initial transient growth, they are rolled up near the centre of the pipe, approaching the unstable eigenfunctions. The exponential growth rate is in good agreement with the linear stability analysis.
Figure 12. Isosurfaces of the streamwise velocity of the NLOP at (a) 
$C=4$ and (b) 
$C=5$, at 
$Re=3000$, 
$\mathcal {T}=400$. Red/yellow are 
$30\,\%$ of the min/max value.
For the interesting case at 
$C=6$, an NLOP with structure similar to that found at 
$C<5$ was not found for similar 
$E_0$ and 
$\mathcal {T}$. By reducing the initial energy and extending the target time substantially (
$E_0 \sim O(10^{-7})$, 
$\mathcal {T}=1000$), we finally located a new type of NLOP. (Note that the energy of the convective state itself is also much reduced compared to shear turbulence, seen in figure 4(a).) Energy growth as a function of time is shown in figure 13. Optimals at all energies experience a period of energy growth at the beginning. Then for a small initial energy 
$E_0=1\times 10^{-8}$, it decays exponentially. At larger initial energy, below the critical energy, the energy grows a little but finally still decays, albeit slowly. For this reason, the long optimisation time 
$\mathcal {T}=1000$ was required. At a slightly larger 
$E_0$ still, the critical energy jump is observed; see figure 13(b).
Figure 13. Optimal energy growth for 
$C=6$, 
$Re=3000$, 
$\mathcal {T}=1000$.
Isosurfaces of the optimal at different 
$E_0$ for 
$C=6$ are shown in figure 14. Unlike the NLOP at 
$C=0\unicode{x2013}4$, its structure is observed to change as 
$E_0$ is increased. The linear optimal is shown in figure 14(a), which is distributed regularly in the streamwise direction. It is very different from the linear optimal of the isothermal case, which consists of streamwise rolls, and is different from the convective eigenfunctions of figure 8. For larger 
$E_0$, the NLOP shown in figure 14(b) has a more spiral structure. As the initial energy increases further in figures 14(c,d), it localises in the streamwise direction. The last of these, figure 14(d), is the minimal seed that eventually triggers turbulence. Its cross-section in figure 15 shows two pairs of very weak streamwise vortices in the centre of the pipe, with very small magnitude compared to the streamwise perturbation. The dynamics of this optimal is investigated in the next subsection.
Figure 14. Isosurfaces of streamwise perturbation at 
$C=6$. Linear optimal (a) 
$E_0=1\times 10^{-8}$, and NLOPs at (b) 
$E_0=1\times 10^{-7}$, (c) 
$E_0=3\times 10^{-7}$, (d) 
$E_0=4\times 10^{-7}$. Red/yellow are 30 % of the min/max.
Figure 15. Cross-section of NLOP at 
$C=6$, 
$E_0=4\times 10^{-7}$. Contours for streamwise perturbations, while the arrows represent the small cross-stream components; the largest arrow has magnitude 
$4.44\times 10^{-7}$.
3.4. Transition to convective turbulence
Transition to shear turbulence via the NLOP has been widely researched for isothermal flow (Pringle & Kerswell Reference Pringle and Kerswell2010; Pringle et al. Reference Pringle, Willis and Kerswell2012; Marensi et al. Reference Marensi, Willis and Kerswell2019), therefore we focus on the transition to convective turbulence (
$C\geqslant 4$) here. In figure 16(a), transition at 
$C=4$ and 
$Re=3000$ is compared for several initial conditions – shear turbulence, the minimal seed, and unstable eigenfunctions. Here, 
$E_{3d}$ is the energy of the streamwise-dependent component of the flow (modes 
$k\ne 0$ in the Fourier expansion (2.23)). For the minimal seed, 
$E_{3d}$ increases substantially at first, reaching the shear turbulence state. It then decays to the convective turbulence state. This implies that shear turbulence can still be triggered in the buoyancy regime, but it is not sustained. Starting with a shear turbulence state taken from an isothermal simulation, there is an immediate decay at first, then 
$E_{3d}$ grows exponentially towards convective turbulence. Interestingly, starting from the unstable eigenfunctions, the dynamics approaches a travelling wave solution (mode 2) or relative periodic solutions (modes 1 and 3). These states can be approached for a long time, suggesting that the travelling wave and relative periodic orbit solutions are only weakly unstable. An equilibrium state and periodic motion have also been reported in early research (Kemeny & Somers Reference Kemeny and Somers1962; Scheele & Hanratty Reference Scheele and Hanratty1962; Yao Reference Yao1987), but at low Reynolds number (
$Re<2000$).
Figure 16. Plots of 
$E_{3d}$ (energy of the streamwise-dependent component of the flow) for different initial conditions at (a) 
$C=4$ and (b) 
$C=5$, at 
$Re=3000$. There are three unstable modes at 
$C=4$: mode 1 (
$m=1$, 
$k=2$), mode 2 (
$m=1$, 
$k=1$) and mode 3 (
$m=1$, 
$k=3$). There is only one unstable mode at 
$C=5$ (
$m=1$, 
$k=1$).
At 
$C=5$ (figure 16b), the unstable eigenfunction does not lead to a travelling wave or periodic solution, and instead transitions to convective turbulence directly. This may be due to the enhanced instability of these equilibrium solutions at larger 
$C$. The minimal seed at 
$C=5$ also triggers the convective turbulence directly, or at least, there is no period of decay after the initial growth, which would indicate a distinct switch from shear turbulence to convective turbulence. Meanwhile, starting from shear turbulence, there is a clear initial decay, followed by exponential growth at a rate that matches the value starting from the eigenfunction.
The transition process at 
$C=6$, where the laminar state is linearly stable, is quite different, shown in figure 17. Compared with transition via the NLOP in the isothermal case (Pringle et al. Reference Pringle, Willis and Kerswell2012), the transition process is much slower: the earliest time turbulence is seen at approximately 
$t=1000$ in figure 17, hence the large target time 
$\mathcal {T}=1000$ required. The edge state (Duguet, Willis & Kerswell Reference Duguet, Willis and Kerswell2008), here to the convective turbulence state, appears to be less unstable than that for the shear turbulence state, as intermediate energies before transition can be achieved for long times with little refinement of the initial energy. Some typical coherent structures that appear in the process of transition are presented in figure 18, starting from slightly different initial energies, 
$E_0=4\times 10^{-7}$ (left-hand plots) and 
$E_0=4.2\times 10^{-7}$ (right-hand plots). Initially, the contours of minimal seed are tightly layered and forward-facing, but because the base profile is M-shaped, the structures are inclined into the shear; see figure 14(d). (For the isothermal minimal seed, similar layers are inclined into the shear, but without the region of reversed flow, they are backward-facing and located closer to the wall (Pringle et al. Reference Pringle, Willis and Kerswell2012).) Therefore, between 
$t=0$ and 
$t=T_1$, there is a great energy growth in a short time through the Orr mechanism (Orr Reference Orr1907; Pringle et al. Reference Pringle, Willis and Kerswell2012). The flow then evolves to a simpler organised state at 
$t=T_2$ (figure 18b), and the flows still look similar. At 
$t=T_3$ (figure 18, right-hand plot), the structure is more elongated, and looks similar to the travelling wave solution of figure 19(a). For slightly lower initial energy, the state gradually decays exponentially to 
$T_4$ (figure 18d, left-hand plot) with a very small decay rate. For the slightly larger initial energy, it evolves into the convective state at 
$T_4$ (figure 18d, right-hand plot). A relatively rapid increase of energy is observed at approximately 
$T_4$ for this case (
$E_0=4.2\times 10^{-7}$), and a similar final growth stage is observed for the NLOPs of larger initial energies, at earlier times, seen in figure 17. Overall, the transition process of the minimal seed is complex and long compared to that of the minimal seed for transition to shear turbulence.
Figure 17. Time series of energy starting from the NLOP obtained at several 
$E_0$ values, at 
$C=6$, 
$Re=3000$. The cases 
$E_0=4\times 10^{-7}, 4.6\times 10^{-7}, 5\times 10^{-7}$ were also shown in figure 13(a) up to 
$t=\mathcal {T}=1000$.
Figure 18. Isosurfaces of streamwise velocity starting from NLOPs for 
$E_0=4\times 10^{-7}$ in the left-hand plots, and 
$E_0=4.2\times 10^{-7}$ in the right-hand plots, at times (a) 
$t=T_1$, (b) 
$t=T_2$, (c) 
$t=T_3$, (d) 
$t=T_4$, for 
$C=6$, 
$Re=3000$. Times are marked in figure 17. Red/yellow are 30 % of the min/max value.
Figure 19. Travelling wave solution approached during transition at 
$C=4$, 
$Re=3000$. Isosurfaces of (a) streamwise perturbation and (b) streamwise vorticity. Red/yellow are 30 % of the min/max value.
3.5. Travelling wave and periodic orbit solutions
In recent years, a series of invariant solutions of the Navier–Stokes equations has been found for pipe flow (Hof et al. Reference Hof, Van Doorne, Westerweel, Nieuwstadt, Faisst, Eckhardt, Wedin, Kerswell and Waleffe2004; Wedin & Kerswell Reference Wedin and Kerswell2004; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007) and other shear flows (Nagata Reference Nagata1990; Toh & Itano Reference Toh and Itano2003; Reetz, Kreilos & Schneider Reference Reetz, Kreilos and Schneider2019). These invariant solutions have significantly enriched our understanding of the transition and maintenance of turbulence. A subset of these invariant solutions is found in the edge between the laminar and turbulent attractors, and sheds some light on the transition to turbulence (Duguet et al. Reference Duguet, Willis and Kerswell2008; Avila et al. Reference Avila, Barkley and Hof2023). Within turbulence, a well-known periodic solution for plane Couette flow (Kawahara & Kida Reference Kawahara and Kida2001) successfully reproduces the self-sustaining process of shear turbulence (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995). Although these solutions are unstable, the dimensions of their unstable manifolds in phase space are typically low (Kawahara Reference Kawahara2005; Kerswell & Tutty Reference Kerswell and Tutty2007). A generic turbulent state may approach them and spend a substantial fraction of its lifetime in their neighbourhood (Kawahara, Uhlmann & Van Veen Reference Kawahara, Uhlmann and Van Veen2012).
In figure 16, the dynamics approached a travelling wave solution and periodic solutions. Using states from the trajectories seen in figure 16(a) at later times as initial guesses, the Jacobian-free Newton–Krylov method (Willis Reference Willis2019) converged very quickly to a travelling wave solution and periodic solutions. The travelling wave solution is plotted in figure 19, which has wave speed 
$c=0.6783$. Streamwise vorticity is localised in the streamwise direction, located at the bends of the high-speed streak. An initial state was also taken from 
$t=1500$ for 
$E_0=4.2\times 10^{-7}$ of figure 17, and rapidly converged to a periodic orbit with weak variation and period 
$68.4$. Linear instability analysis of this periodic orbit solution shows that it has only one real eigenvalue, with small real part 
${\rm Re}(\sigma )=8.06\times 10^{-4}$, so that excessive refinement of the initial energy is not required to stay close to it for significant time before diverging to the convective turbulence or the laminar state (stable at 
$C=6$). As this periodic orbit has only one unstable direction, it is an attractor within the boundary between the laminar and convective states, i.e. it is an ‘edge state’ (Duguet et al. Reference Duguet, Willis and Kerswell2008).
The periodic solutions approached from the first and third unstable eigenfunctions at 
$C=4$, denoted 
$PO_1$ and 
$PO_2$, respectively, have long periods, 
$T_1=724.4$ and 
$T_2=623.72$. They are similar, only with a different dominant axial wavenumber (
$k=2$ for 
$PO_1$, and 
$k=3$ for 
$PO_2$). We here pay attention to the periodic solution triggered by the first unstable mode. Time series of energies over approximately two periods of 
$PO_1$ are shown in figure 20(a): the energy of the rolls, 
$E_{roll}=E(u_r)+E(u_\phi )$, streaks, 
$E_{streak}=E(u_z)_{k=0}$, and waves, 
$E_{waves}=E(u_z)_{k\ne 0}$, where the subscripts 
$k=0$ and 
$k\ne 0$ indicate streamwise independent and dependent components extracted via the Fourier decomposition. As a single streak forms at the centre of the pipe, we have not split azimuthal dependence through 
$m$. For the well-known periodic process studied in channel flow by Hamilton et al. (Reference Hamilton, Kim and Waleffe1995), the peak of the roll energy closely corresponds to the valley of the streak energy, but for 
$PO_1$, these energies are more offset. Contours of streamwise perturbation velocity 
$u_z$ at 
$t=500,700,900,1100$ are plotted in figure 21, which shows how the structure changes over one period. At 
$t=500$, the streak energy is at its lowest, and the velocity field resembles the unstable eigenfunction, only with a weak break of shift and reflection symmetry. Then by 
$t=700$, the unstable eigenfunction has grown until regions of positive velocity have merged, forming a high-speed streak in the centre of the pipe. At this stage, the wavy high-speed streak in the centre of the pipe is similar to that of the travelling wave of figure 19. In the presence of the streak, the roll perturbations decay at 
$t=900$, then return once the streak has weakened at 
$t=1100$. Therefore, three typical stages are observed in the periodic solution, i.e. growth of the eigenfunction, formation of the streak, and decay of the streak. These processes can be understood from the changes in the mean velocity profile, seen in figure 20(b). At first, when the perturbation is weak, the velocity profile is close to the laminar solution (
$t=500$). Linear instability takes charge of the dynamics, and the unstable eigenfunction grows at this stage. Then, when the perturbation excites a sufficiently strong streak near the axis, the velocity profile is flattened (
$t=700,900$) and the linear instability is suppressed, causing the perturbation to decay. Once the perturbation decays to a certain level, the velocity profile becomes unstable again (
$t=1100$) and the process repeats. (The full periodic process is shown in the supplementary movie available at https://doi.org/10.1017/jfm.2024.589.)
Figure 20. (a) Time evolution of the periodic solution 
$PO_1$ at 
$C=4$, 
$Re=3000$, where 
$E_{roll}=E(u_r)+E(u_\phi )$, 
$E_{wave}=E(u_z)_{k\ne 0}$ and 
$E_{streak}=E(u_z)_{k=0}$. (b) The turbulent mean velocity profile at several times.
Figure 21. Contours of periodic solution 
$PO_1$ at (a) 
$t=500$, (b) 
$t=700$, (c) 
$t=900$, (d) 
$t=1100$, coloured by streamwise velocity disturbance.
In the chaotic convective state at 
$C=4$, phenomena similar to those of the periodic orbits occur, but the laminar flow is unstable to multiple modes. Contours of convective turbulence show that the flow still exhibits the three typical processes, i.e. the growth of unstable modes (figure 22a), the formation of the high-speed streak (figure 22b), and the decay of the streak and reappearance of the unstable modes (figures 22c,d). The scales of the observed modes are in good agreement with the first unstable mode and the third unstable mode (
$k=2$ and 
$k=3$), which is consistent with the flow wandering between the two periodic solutions. The scale of the travelling wave (
$k=1$) is hardly observed.
Figure 22. Contours of transient convective turbulence at (a) 
$t=2500$, (b) 
$t=2750$, (c) 
$t=3000$, (d) 
$t=3250$, picked from figure 16 when the turbulence is triggered by shear turbulence at 
$C=4$. The contours are coloured by streamwise velocity disturbance.
At larger 
$C$, the periodic solutions become more unstable and the flow is more chaotic. However, the typical process can still be observed, as in figure 23, which shows two transient flow states taken from the flow at 
$C=5$. Figure 23(b) shows the typical growth of an unstable mode upon a mean velocity profile of figure 23(a). It is observed that the profile is close to the laminar profile, but does not go so close as the case at 
$C=4$ of figure 20(b), as the laminar solution is also more unstable. Figure 23(c) shows the high-speed streak formed by growth of the unstable mode, whose profile is almost flattened and loses linear instability. These typical flow states can often be observed at 
$C=5$, and can even be observed at 
$C=10$, although less clearly. Therefore, it appears that the convective turbulence is self-sustained by these typical stages, i.e. the growth of unstable modes, the formation of streaks and the decay of streaks.
Figure 23. (a) The turbulent mean profiles of the flow at (b,c). Contours of transient convective turbulence at (b) 
$t=1150$, (c) 
$t=1350$, picked from figure 16 when the turbulence is triggered by shear turbulence at 
$C=5$. The contours are coloured by streamwise velocity disturbance.
4. Conclusion
In this study, a new direct numerical simulations model that includes time-dependent heat flux and background temperature gradient is established. The results show good consistency with the experiment, and a little improvement over the model of Marensi et al. (Reference Marensi, He and Willis2021), which also had time-dependent heat flux, but a spatially uniform heat sink term. Simulations at different values of the buoyancy parameter 
$C$ confirm three typical flow states, i.e. shear turbulence, laminarisation and convective turbulence. The distribution of these states in the parameter space is consistent with the calculations of Marensi et al. (Reference Marensi, He and Willis2021). Detailed examination of the convective turbulence found that it lacks near-wall rolls, which enhance the mixing of fluid near the wall and centre, therefore the heat transfer is seriously reduced. The linear stability analysis further verifies that the convective turbulence is driven by linear instability (Yao Reference Yao1987; Su & Chung Reference Su and Chung2000; Khandelwal & Bera Reference Khandelwal and Bera2015; Marensi et al. Reference Marensi, He and Willis2021). The structures of unstable eigenfunctions are found to be concentrated mainly on the centre of the pipe, which explains the similar flow characteristics observed in the transient convective turbulence.
Nonlinear non-modal stability analysis is extended to the heated pipe case to determine the effects of buoyancy on the smallest perturbation that triggers transition, the ‘minimal seed’. As 
$C$ is increased from 
$0$, it is found that the structure of the minimal seed becomes thinner, rolls move nearer to the wall, and the critical energy for transition increases. Maximum energy growth occurs earlier and is greatly reduced. Importantly, this suggests that strategies for exciting shear turbulence to encourage greater heat flux do not carry over straightforwardly from the isothermal to the heated pipe case. The branch of the NLOP that triggers turbulence in isothermal flow could not be tracked beyond 
$C=5$, due to the combined reduction in growth and enhancement of linear instability. Instead, a new type of NLOP arises for a longer target time and a small initial energy. The new NLOP changes structure as the initial energy is increased, localising in the streamwise direction. It is concentrated towards the centre of the pipe, triggering convective turbulence directly, but over a long period.
The transitions towards convective turbulence at 
$C=4,5,6$ triggered by the NLOP were compared with transition via the unstable eigenfunctions. At 
$C=4$, it is found that starting from the unstable eigenfunctions leads to a travelling wave solution or periodic solutions, while the minimal seed first triggers shear turbulence, then the shear turbulence decays to convective turbulence. At 
$C=5$, both unstable eigenfunctions and the minimal seed trigger convective turbulence directly. At 
$C=6$, the transition starting from the new NLOP is much slower. The edge state to the convective turbulence appears to be less unstable than that of shear turbulence, and intermediate energies before transition can be achieved for a long time with little refinement of the initial energy. The approached travelling wave solution and periodic solutions were calculated using the Jacobian-free Newton–Krylov method, and their stability calculated using Arnoldi iteration. The periodic orbit solution at 
$C=6$ has only one unstable direction, as expected for an edge state.
The periodic solution is found to distil three typical processes, i.e. the growth of an unstable eigenfunction, the formation of streaks and the decay of the streaks. Analysis of the mean velocity profile showed that the periodic process is caused by the appearance and suppression of linear instability of the mean velocity profile. This is fundamentally different from the self-sustaining process of isothermal shear flow, which occurs in the absence of linear instability of the mean flow. Chaotic convective turbulence at 
$C=4$ is still dominated by the three typical processes, but consists of the scales of two periodic solutions, consistent with the dynamics wandering between the two solutions in phase space. A similar phenomenon is still often observed at 
$C=5$, and can be recognised even up to 
$C=10$.
We caution that of our work has focused on the case 
$Re=3000$. However, the form of the minimal seed is known to be robust in the isothermal case with respect to large changes in 
$Re$ and flattening of the base profile (Marensi et al. Reference Marensi, Willis and Kerswell2019). Due to the collapse with respect to the buoyancy parameter 
$Bo$ for the transition seen in figure 3(b), we expect similar behaviour to be exhibited at larger 
$Re$, in particular, for values of 
$C$ around where transition from shear turbulence to the convective state is observed.
Supplementary movie
A supplementary movie is available at https://doi.org/10.1017/jfm.2024.589.
Acknowledgements
This work used the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1).
Funding
S.C. acknowledges funding from Sheffield–China Scholarships Council PhD Scholarship Programme (CSC no. 202106260029).
Declaration of interests
The authors report no conflict of interest.
