Suppression of turbulence and travelling waves in a vertical heated pipe

Turbulence in the flow of fluid through a pipe can be suppressed by buoyancy forces. As the suppression of turbulence leads to severe heat transfer deterioration, this is an important and undesirable phenomenon in both heating and cooling applications. Vertical flow is often considered, as the axial buoyancy force can help drive the flow. With heating measured by the buoyancy parameter $C$, our DNS show that shear-driven turbulence may either be completely laminarised or transitions to a relatively quiescent convection-driven state. Buoyancy forces cause a flattening of the base flow profile, which in isothermal pipe flow has recently been linked to complete suppression of turbulence (K\"{u}hnen et al. Nat. Phys., 2018), and the flattened laminar base profile has enhanced nonlinear stability (Marensi et al. JFM, 2019). In agreement with these findings, the nonlinear lower-branch travelling-wave solution analysed here, which is believed to mediate transition to turbulence in isothermal pipe flow, is shown to be suppressed by buoyancy. A linear instability of the laminar base flow is responsible for the appearance of the relatively quiescent convection driven state for $C\gtrsim 4$ across the range of Reynolds numbers considered. In the suppression of turbulence, however, i.e. in the transition from turbulence, we find clearer association with the analysis of He et al. (JFM, 2016) than with the above dynamical systems approach, which describes better the transition to turbulence. The laminarisation criterion He et al. propose, based on an apparent Reynolds number of the flow as measured by its driving pressure gradient, is found to capture the critical $C=C_{cr}(Re)$ above which the flow will be laminarised or switch to the convection-driven type. Our analysis suggests that it is the weakened rolls, rather than the streaks, which appear to be critical for laminarisation.


Introduction
Most energy systems rely on fluids to transfer heat from one device to another to facilitate power generation, provision of heating or production of chemicals. Flows are often forced through channels or arrays of pipes taking heat away from the surfaces. In a nuclear reactor, for example, the reactions occur within the fuel pins, which are cooled by flow of coolant through the channels formed by arrays of fuel pins to maintain their temperature within a specific limit as well as transferring energy to the steam generators. In an isothermal flow, the volume flux is driven by an externally applied pressure gradient, and the flow is referred to as 'forced'. In a vertical configuration, however, buoyancy resulting from the lightening of the fluid close to the heated wall can provide a force that partially or fully drives the flow, referred to as mixed or natural convection, respectively. When heat flux is very high, we can have a 'supernatural' state of flow, where the buoyancy is sufficiently strong that a reversed pressure gradient may be necessary to limit or maintain a constant volume flux. Under certain conditions (e.g. the Boussinesq approximation) an upward heated flow may be considered equivalent to a downward flow cooled at the boundary (Appendix A).
Mixed convection is of significant importance to engineering design and safety considerations and as such extensive research has been carried out to develop engineering correlations (Jackson et al. 1989;Yoo 2013), turbulence models (Kim et al. 2008;Bae 2016) and a better understanding of the physical flows (You et al. 2003). A particularly interesting physics is that the flow, at a Reynolds number where shear-driven turbulence is ordinarily observed, in the presence of buoyancy may be partially or fully laminarised, or becomes a convection-driven turbulent flow (i.e. natural convection, referred to above). Heat transfer may be significantly impaired under such conditions. He et al. (2016) (hereinafter referred to as HHS) modelled the effect of buoyancy using a prescribed body force, with linear or step radial dependence, without solving the energy equation. They attributed the suppression of turbulence to a reduction in the apparent Reynolds number of the flow, as measured by the pressure gradient required to drive the flow. Thus, the forced flow is compared with the unforced "equivalent pressure gradient" reference flow.
Meanwhile, in ordinary (isothermal) pipe flow, Kühnen et al. (2018), observed relaminarisation attributed to flattening of the base flow profile. The idea of flattening was first suggested by Hof et al. (2010) who showed that when two puffs were triggered too close to each other the downstream puff would collapse due to the flattened streamwise velocity profile induced by the upstream puff. In the experiments of Kühnen et al. (2018) the flattening was introduced by a range of internal and boundary flow manipulations and a full collapse of turbulence was obtained for Reynolds numbers up to 40 000. Marensi et al. (2019) showed the complement effect, i.e. the enhanced nonlinear stability of the laminar flow. They found that the minimal seed (smallest amplitude disturbance) for transition is 'pushed out' from the laminar state to larger amplitudes when the base flow is flattened, thus implying that a flattened base profile is more stable than the parabolic profile. Here, buoyancy forces also have a flattening effect and turbulence may be partially or fully suppressed. Furthermore, early experimental observations (Hanratty et al. 1958;Kemeny & Somers 1962;Scheele & Hanratty 1962) and subsequent linear (Yao 1987a,b;Yao & Rogers 1989;Chen & Chung 1996;Su & Chung 2000) and weaklynonlinear (Rogers & Yao 1993;Khandelwal & Bera 2015) stability analyses suggested that, for sufficiently large heating, the flow becomes unstable and transitions to a new non-isothermal equilibrium state characterised by large-scale regular motions. In agreement with the experiments, the linear theory showed that this instability can occur at low Reynolds number (below 100) and for Re > 50 the critical value of the Rayleigh number is almost independent of Re (Yao 1987a). The first azimuthal mode was found to be the least stable (Yao 1987a;Su & Chung 2000), consistent with the double-spiral patterns observed experimentally (Hanratty et al. 1958) and the instability was linked to the inflectional velocity profile in the buoyancy-assisted case. As suggested by Su & Chung (2000), a competition between different mechanisms -driven by either shear or convection -thus exists and understanding its effect on the nature of the flow is the object of our study.
In particular, in this work, we are interested in whether a flow is turbulent or laminar under certain heating conditions and when a turbulent flow may be laminarised or vice versa under the influence of buoyancy. We address this question for a vertically heated pipe, initially in the dynamical systems context through linear stability and by investigating how travelling wave solutions are affected by the buoyancy force. Next, the nature of the laminarisation is considered. In isothermal flow at transitional Reynolds numbers, the shear-driven state is known to be metastable -the probability of laminarisation follows a Poisson process with a laminarisation rate that depends on the Reynolds number.
In any practical setting, where a pipe is of finite length, its length affects the probability of turbulence surviving to the end of a pipe. Hence a range of Reynolds numbers for transition are quoted, typically between 2000 and 2300. Therefore, we do not attempt to quantify the full statistical nature of the transition in the heated case, but instead we focus on the phenomenological-based 'equivalent pressure-gradient' analysis of HHS. Through the above approaches, i.e. linear stability, nonlinear travelling-wave and 'equivalent pressure-gradient' analyses, we aim to elucidate the physical mechanisms underlying the buoyancy-suppression of turbulence, illustrating the bistability nature of such flows.

Nonlinear dynamical systems view
In subcritical wall-bounded shear flows, turbulence arises despite the linear stability of the laminar state (Drazin & Reid 2004;Schmid & Henningson 2001). The implication is that the observed transition scenario can only be triggered by finite amplitude disturbances. In the last 30 years our understanding of transition to turbulence in such flows has greatly benefited from a fully nonlinear geometrical approach which adopts concepts from the dynamical systems theory. In this view, the flow is analysed as a huge (formally infinite-dimensional) dynamical system in which the flow state evolves along a trajectory in a phase space populated by various invariant solutions, travelling waves (TWs) and periodic orbits (POs). Nonlinear travelling wave solutions were first discovered numerically in the early 1990s for plane Couette flows (Nagata 1990) and in the 2000s for pipe flows (Faisst & Eckhardt 2003;Pringle & Kerswell 2007). Since then, partly thanks to the advances in our computational and experimental capabilities, a growing amount of evidence has been collected for their dynamical importance (see reviews Kerswell 2005;Eckhardt et al. 2007;Kawahara et al. 2012;Graham & Floryan 2021). These solutions, often referred to as "exact coherent states/structures" (ECSs), are believed to act as organising centres (Waleffe 2001) in phase space, in the sense that, when the flow state approaches them, spatio-temporally organised patterns (streaks and streamwise rolls) are observed (Hof et al. 2004;Kerswell & Tutty 2007). ECSs are finite-amplitude non-trivial solutions disconnected from the laminar state and enter via saddle-node bifurcations as the flow rate is increased. Some solutions, typically those of higher spatial symmetry, exist at flow rates much below that at which transition is usually observed (Pringle et al. 2009). ECSs are linearly unstable, although with only a few unstable directions. They may be divided into 'upper-branch' and 'lower-branch' states, depending on whether they are associated with a high or low friction factor. Lower branch solutions are representative of the laminar-turbulent boundary -the so called "edge of chaos" (Itano & Toh 2001;Schneider & Eckhardt 2006) -which separates initial conditions that lead to turbulence from those that decay and relaminarise. The edge comes closest to the laminar equilibrium at the "minimal seed" for transition (Kerswell 2018). Lower-branch solutions are believed to mediate the transition to turbulence (Duguet et al. 2008;Schneider et al. 2007), while some upper-branch solutions are embedded in the turbulent attractor and are representative of the turbulent dynamics .
Here, we are interested in studying how travelling wave solutions are affected by the buoyancy force in a vertical heated pipe, and, in analysing their dynamics, we aim to elucidate the physical mechanisms underlying the buoyancy-suppression of turbulence. The transition between regimes is first investigated using linear stability in §3.2, followed by analysis of travelling waves in §3.3.

Equivalent pressure-gradient (EPG) analysis of HHS
Rather than simulating a temperature field, to reduce complexity HHS considered a fixed radially-dependent axial body force that models the buoyancy force, and applied this to isothermal flow. Conventionally, heated flows are compared with the isothermal (unforced) flow at equivalent flow rate (EFR), but HHS observed better comparison with flows at the equivalent pressure gradient (EPG). In particular, after careful analysis, they observed that adding the radially-dependent force does not alter the turbulent viscosity of an unforced flow driven by the same pressure gradient (see figure 10 therein). The unforced EPG flow is therefore a reference flow for cases with the extra radiallydependent forcing. Note that in a fixed mass-flux calculation, the pressure gradient reduces in response to driving from the buoyancy. Given a heated flow at a particular Reynolds number Re (defined in terms of the mass flux), to determine the Reynolds number of the EPG flow, one must split the mass flux into contributions from the pressure gradient and from the buoyancy. The former component determines the 'apparent Reynolds number' Re app of the EPG flow. Laminarisation of the body forced flow is observed to occur when its Re app is consistent with the Re at which laminarisation occurs in isothermal flow. Further details of the analysis are provided in §3.4 and HHS prediction is compared with a suite of simulations in §3.5.

Formulation
Consider a vertically aligned circular pipe of diameter D, with the flow of fluid upwards. We model a short pipe section of length L (figure 1(left)) and let {u(x, t), p(x, t), T (x, t)} be the velocity, pressure and temperature fields, respectively. The fluid has kinematic viscosity ν, density ρ, volume expansion coefficient γ and thermal diffusivity κ. Under the Boussinesq approximation, density variations are ignored except where they appear in terms multiplied by the acceleration due to gravity, gẑ, leading to the governing equations where T ref is a reference temperature defined in the following subsection, and d z P is the pressure gradient for laminar flow with bulk velocity U b . We suppose that U b is fixed, in which case β(u) adjusts to maintain fixed bulk velocity. We also suppose that the temperature of the wall T w and the bulk temperature T b are fixed. The latter is achieved by including a uniform heat sink (t) which adjusts to maintain the fixed bulk value T b . For such a flow, we can introduce axial periodicity, so that (t) may be considered equivalent to the rate at which heat absorbed by the fluid would otherwise be carried out of the section of pipe †. For laminar flow, the flow is purely axial so that radial heat transport is purely conductive. If 0 is the heating rate for the laminar case, then the observed quantity Nu :=¯ / 0 is the Nusselt Number, where the overbar (•) denotes time average.

Non-dimensionalisation
Given the temperature at the wall T w and the bulk temperature T b , we put ∆T = 2(T w − T b ) and take a reference temperature where T c is the centreline temperature for the case of laminar flow. (The choice for T ref does not influence the flow, since the constant γ g T ref could be absorbed into the pressure gradient.) We introduce the dimensionless temperature Θ = (T − T c )/∆T . Let the pipe radius R = D/2 be the length scale and the isothermal laminar centreline velocity 2 U b be the velocity scale. The corresponding time scale is thus R/(2 U b ). Hereafter, all variables are dimensionless except (t) which always appears in the dimensionless ratio / 0 , i.e. the instantaneous Nusselt number. Non-dimensionalising with these scales, for the temperature equation we find For the laminar case, Θ = Θ lam = r 2 , we find (2.5) † Spatial periodicity limits the domain over which wall friction is averaged, which can lead to unrealistic fluctuations (mean-square variations from the time average) in the bulk velocity. We therefore assume constant flux.
Plugging this ∆T back in to (2.4), we obtain the dimensionless temperature equation where Re := U b D/ν is the Reynolds number and P r := ν/κ is the Prandtl number. For the momentum equation we find The coefficient of the buoyancy term can be written Although the Grashof number is in common use, from Gr it is difficult to judge the magnitude of the buoyancy force relative to the pressure gradient of the laminar flow for this particular configuration. We therefore write the dimensionless momentum equation as where C measures the buoyancy force relative to the force that drives the laminar isothermal shear flow, C = Gr/(4 Re 2 ) 4/Re := Gr 16 Re . (2.10) The laminar velocity and laminar temperature profiles for this configuration are (2.11) and the no-slip and fixed-temperature boundary conditions at r = 1 are u = 0, Θ = 1, (2.12) respectively, while periodic boundary conditions are applied in the streamwise direction. The laminar velocity profiles for different C are shown in figure 1(right). The isothermal pipe flow is recovered for C = 0 (no buoyancy force) and P r = 0 (temperature diffuses immediately), with the parabolic laminar profile U 0 = 1 − r 2 . For a statistically steady flow, Reynolds averaging is both time averaging and cylindrical surface averaging, where the latter is denoted as (2.13) Turbulent fluctuations are calculated as deviations from the mean components of the

Numerics
Simulations were carried out using the Openpipeflow solver (Willis 2017), modified to include timestepping of the temperature field and the buoyancy term in the momentum equation. A variable q(r, θ, z) is discretised using a non-uniform grid in the radial direction with points clustered near the wall and Fourier decompositions in the azimuthal and streamwise directions, namely where α = 2π/L is the streamwise wavenumber and m p determines the azimuthal periodicity (m p = 1 for no discrete rotational symmetry). Radial derivatives are evaluated using central finite differences with a nine-point stencil. At Re = 5300, in a L = 5D long pipe we use a spatial resolution of (N × M × K) = (64 × 96 × 96), which ensures a drop of at least 4 orders of magnitude in the spectra and provides the correct value for the friction factor, as reported in the literature (Eggels et al. 1994). Following the 3/2 dealiasing rule, variables are evaluated on an N × 3M × 3K grid in physical space. A second-order predictor-corrector scheme is employed for temporal discretisation, and a fixed timestep of 0.01 is used. This is sufficient to ensure that the time discretisation error is no larger than the spatial discretisation error (measured by the corrector and spectra respectively) and corresponds to a CFL-number of approximately 0.2 − 0.25. Data for simulations for various Gr = 16 Re C and constant Re = 5300, P r = 0.7 are shown in figure 2. There is good agreement with numerical data (You et al. 2003) and experimental data (Steiner 1971;Carr et al. 1973;Parlatan et al. 1996).

Travelling wave solutions
In order to apply dynamical systems theory, the discretised momentum and temperature equations are formulated as an autonomous dynamical system (Viswanath 2007;Willis et al. 2013): where X is the vector of dependent variables, here X = (u, Θ), and p is the vector of parameters of the system, p = (Re, C). The simplest solution is an equilibrium, which satisfies F(X; p) = 0. For pipe flow, the only equilibrium solution is the laminar solution. Travelling wave solutions satisfy X(x, t) = g(ct) X(x, 0), where here g(l) applies a streamwise shift by l, and c is the phase speed. Travelling waves are also known as 'relative' equilibrium solutions, as they are steady in a co-moving frame. They therefore satisfy for some vector (X, l, T ), and hence can be calculated via a root solving method. The most popular method at present is the Newton-Krylov method. (Note that in addition to (2.16), two extra constraints are required to match the extra unknowns l, T ; see Viswanath (2007).) Time-dependent periodic orbits may also be calculated by this method. Typically periodic orbits originate via a Hopf bifurcation off a travelling wave, but are not discussed further in this work. Stability of the solutions is calculated using the Arnoldi method to solve the eigenvalue problem where σ is the growth rate and the operator on the right hand side is linearised about the travelling wave X 0 by taking ||dX|| The Newton-Krylov and Arnoldi solver, already available as a utility of Openpipeflow (Willis 2017), were integrated with the time-stepping code described in §2.2 for heated pipe flow.

Results and discussion
All results presented herein pertain to the case P r = 0.7 and constant volume flux. This relatively low Prandtl number is a reasonable starting choice for the applications we are interested in, where most gasses have P r ≈ 0.7, e.g. CO 2 . In large scale cooling applications using liquid metal, P r is much smaller. Cases where P r > 1 (e.g. P r = 7 for water) are more expensive numerically due to a need for higher resolution for the temperature field.

Direct Numerical Simulations
Simulations were performed in a pipe of length L = 5D for a range of Reynolds numbers to study the effect of the buoyancy parameter C. Results are first shown for a relatively low Reynolds number, Re = 2500. Figure 3 shows complete relaminarisation of transitional turbulence in response to the introduction of buoyancy for intermediate values of . Relaminarisation events are revealed by monitoring the energy of the streamwise-dependent component of the flow, denoted E 3D , which shows a rapid decay when the flow relaminarises, E 3D → 0 and ε(t)/ε 0 → 1. At larger C O(10), turbulent fluctuations are not completely suppressed. Instead a convection-driven flow is set up, which becomes stronger as C is increased.
At Re = 5300 the effect of buoyancy is found to be slightly different -turbulence is not observed to undergo complete relaminarisation, but instead transitions directly to a weak convection-driven state. Figure 4 shows simulations with C = O(1) − O(10). The buoyancy causes suppression of the turbulence and therefore a drop in ε(t)/ε 0 , so that the Nusselt number Nu =ε/ε 0 reduces substantially. The corresponding velocity and temperature mean profiles, u z (r) and Θ (r), are shown in bottom graphs of figure 4 together with the laminar profiles at C = 0 for comparison. Cases where turbulence is suppressed exhibit a flattened base velocity profile.
The case for C = 7.5 is shown for longer time in figure 5(left). The shear-driven turbulent state is metastable only, and around t ≈ 2000 turbulence is more suppressed as  there is a switch to the more quiescent convection-driven state. As C is increased further the buoyancy starts to drive a more turbulent convection-driven state. For these cases the velocity profile is more 'M-shaped' as seen in figure 5(right).
The convective state at Re = 5300 and C = 7.5 is visualised in figure 6 together with the metastable shear-driven turbulent state. When comparing the deviations from the are omitted for all trajectories and the curves corresponding to C 12.5 are shifted in time by an arbitrary offset, for clarity only. Right: Snapshots of the mean streamwise velocity profiles u z (r) for the same values of C shown on the left. All the snapshots are taken at t = 1000. For C = 7.5 an additional snapshot (solid grey line with dots) is shown corresponding to t = 7500 (marked with a grey dot on the corresponding trajectory on the left). The thick light-grey line on the right corresponds to the laminar streamwise velocity profile (2.11) with C = 0.
isothermal laminar profile (a-d), both the shear and convective states show a deceleration in the core and acceleration close to the wall, with the convective states showing very smooth and almost z− and θ−independent contour levels. Deviations from the mean profile (e-h), however, reveal that the convective state has larger and more elongated flow structures compared to the shear-driven turbulence. In both types of visualisations it is clear that the small-scale turbulent eddies are strongly suppressed in the convectiondriven flow. Figure 7 shows the type of state seen in simulations, laminar flow (L), shear-driven turbulence (S) and convection-driven flow (C), for a range of Re and C. The initial condition for each simulation was a previously calculated shear-driven state at similar Re. (This is except for Re 2000 and C > 3, where it is clear that the shear-driven state decays immediately, i.e. only the convective state could be supported, and hence the initial condition was of convection type). For each simulation it is relatively easy to distinguish between the shear-and convection-type flows, since the former shows far more chaotic time series and higher heat flux. The case for C = 7.5 in figure 5(left) shows this difference, and also that multiple behaviours are possible at the same parameters for significant periods of time. The shear-driven state is marked if observed for 1000 time units. (It is stable or at least metastable with a long expected lifetime.) A relaminarisation is marked if the energy of the streamwise component of the flow drops below 10 −5 . Overall, figure 7 indicates that as C (or Gr) increases, a larger Re is needed in order to drive shear turbulence, or, equivalently, as Re increases, shear-driven states persist to larger C. For C 4 simulations suggest that a convective instability kicks in, roughly independently of the Reynolds number over this range. In between, it is possible to completely relaminarise flow up to Re ≈ 3500, but at larger Re the progression is as in figure 5 -from a shear-driven turbulent state to a weak convection-driven state, then to a more turbulent convection-driven state as C is increased.
In the following sections we determine whether the boundaries of stability observed

Linear stability analysis
As the transition to shear-driven turbulence in isothermal flow occurs in the absence of a linear instability, this section relates to the transition to convection-driven flow states, in particular with respect to loss of stability of the modified laminar base profile (2.11) for non-zero C. Linear stability of mixed-convection pipe flow has been studied by Yao (1987a); Su & Chung (2000), where the model differs slightly in the boundary condition and form of the heat sink. Our figure 2 suggests these differences make little difference to transition, however, we check for consistency with the nonlinear results of §3.1. As our code uses Fourier expansions in the periodic dimensions, to calculate the eigenfunctions and stability of the base flow (2.11) we need simulate only using a few Fourier modes. The Arnoldi method is employed to accelerate convergence and to access eigenvalues beyond the leading one. Linear stability analysis is performed for azimuthal wavenumbers m = 0, 1, 2 and two streamwise wavenumbers α = 0.628 and α = 1.7 (commensurate with the pipe lengths L = 5D and L = 1.85D used in our DNS study of §3.1 and in the travelling wave analysis of §3.3).
The neutral curves, where the growth rate (σ) = 0, are shown in figure 8. As expected (and as also reported by Yao (1987a); Su & Chung (2000)), the first azimuthal mode is found to be the least stable, it corresponds to the spatially largest mode and is the only mode that can exhibit flow across the axis. (The axisymmetric mode m = 0 is included in the numerical calculations for stability of the m = 1 mode, but we have not observed instability of m = 0 type.) As shown in figure 8, the m = 1 mode exhibits a fairly complex dependence on C, while it is only weakly affected by the axial wavenumber. Indeed, the first branch for α = 0.628 almost coincides with that for α = 1.7 and the other two branches (not shown) are slightly shifted to the right. Consistent with the linear stability of isothermal pipe flow, the critical Reynolds number approaches infinity as C → 0 for any m.
Consistent with the appearance of the convective state found in simulation (figure 7), at C ≈ 4 a linear instability appears, roughly independent of Re for most of the range considered. The corresponding laminar profiles for C = 3−10 are shown in figure 1(right). For C > 4 the profiles present an "M-shape" (independent of Re, see (2.11)), which becomes increasingly more pronounced as C increases. The difference at the centreline is more than 80% for C = 10. The profile at C = 3 is flatter than the parabolic (isothermal) profile, with a centreline difference of almost 30%, but does not have any inflection point. Therefore, in agreement with previous experimental and theoretical studies (Scheele & Hanratty 1962;Yao 1987a;Su & Chung 2000), our analysis suggests that the linear instability of buoyancy-assisted pipe flow is linked to the inflectional velocity profiles ocurring at sufficiently large heating and it is almost independent of Re. Figure 8 also shows that, for C 4, a region of restabilisation is observed as Re is increased. This is also evidenced in figure 9, which shows a region of negative (σ) for 1450 < Re < 6200 at C = 5. Isosurfaces of streamwise vorticity for the eigenfunctions  corresponding to the two neutral points where (σ) becomes positive (Re ≈ 400 and 6200) are also shown in the insets of figure 9. For the larger Reynolds number, Re ≈ 6200, the eigenfunction looks like it is spiralling in the centre and resembles the "spiral" solution found by Senoo et al. (2012), although their visualised solutions are nonlinear.

Continuation from TW N4L
To better understand the effect of buoyancy, we perform a nonlinear analysis, starting from a known TW in isothermal pipe flow (C = Gr = 0) and continuing the solution to larger values. A vast repertoire of TWs has now been compiled in isothermal pipe flows . For our purpose we decided to focus on a fundamental solution, labelled TW N4L (Pringle et al. 2009), which is highly-symmetric (satisfying both shiftreflect and shift-rotate symmetries) and characterised by relatively smooth continuation branches in order to aid the numerical continuation. In Willis et al. (2013), the lower branch of this solution was found to lie on the boundary between the laminar state and turbulence in a 'minimal flow unit'. Localised solutions bifurcate off this class of solutions (Chantry et al. 2014) and are found to mediate transition in extended domains Budanur & Hof 2017).
Following Willis et al. (2016) we start with the 'minimal flow unit' at Reynolds number Re = 2500 with domain (r, θ, z) = [0, 1]×[0, π/2]×[0, 2π/1.7], i.e. m p = 4 and α = 1.7 in (2.14). For isothermal flow (C = Gr = 0), the phase speed of TW N4L is c = 0.61925. The isothermal TW was first reconverged at P r = 0.7 using the Newton solver. A parametric continuation in C to non-zero values was then performed (figure 10) for fixed Re, P r and α. We were able to continue the isothermal solution from C = 0 around positive C and find that it connects with the upper branch at C = 0, then beyond to C ≈ −40. (Negative C corresponds to a downward cooled flow; see Appendix A). As a check, we verified that the values of c = 0.52575 and N u = 2.378 at C = 0 on the upper branch, as well as the mean profiles, matched those of the previously known upper-branch isothermal solution TW N4U with P r = 0.7.
In figure 10(right) it is seen that from C = 0 to C = 6 the Nusselt number Nu increases by approx 0.75. By comparison, along the upper branch, over the large range C = 6 to C = −40, it increases by only a further 1.25. Relatively speaking, the lower branch is rapidly pushed back towards the upper branch over the increase in C and is suppressed altogether for C > 7.5. The mean velocity and temperature profiles at different points along the continuation are shown in figure 11. Observe that the profile in the near-wall region, where rolls and streaks occur, is similar at the saddle-node point (SN) to that  Figure 12(left) shows these rolls (arrows) and streaks (contours) in cross sections of the velocity perturbation at the saddle-node point. The corresponding temperature perturbation field ('thermal streaks') is shown on the right. Similar to its isothermal counterpart, the travelling wave is characterised by fast streaks located near the pipe wall and slow streaks in the interior. The core shows a strongly decelerated region relative to the laminar (isothermal) profile and thus the profile must become steeper at the wall to preserve the mass-flux. The difference from the isothermal TW N4L , however, is less marked in the near-wall region than it is in the core. Continuations were also performed at Re = 2000 and 3000, after reconverging the isothermal TW N4L at these Reynolds numbers. Results are shown in figure 13. The TW survives to larger C as the Reynolds number increases (the saddle-node point of each curve moves to larger C as Re increases). This is consistent with the shear turbulence region in figure 7 persisting to larger C as Re is increased. The saddle-node bifurcations at each Re occur at much larger values of C than those at which suppression of turbulence was observed in the DNS. For example, at Re = 2500 the saddle-node bifurcation occurs at C ≈ 7.5, while in figure 7 shear-turbulence survives only for C 1. This is not so surprising, considering that in isothermal pipe flows the lowest Re at which the N4L travelling-wave solution is found, i.e. Re = 1290 (Pringle et al. 2009), is much below the commonly observed value for transition in experiments (Re ≈ 1800−2300). Furthermore, it should be taken into account that only one TW solution is analysed here -it cannot capture the entire phenomenon of turbulence suppression in a heated pipe flow, although is found to capture some of the fundamental characteristics. Figure 14 shows that, while the lower branch solution for Re = 3000 is on the edge of an attractor for shear-driven turbulence at C = 0, this is no longer the case for C = 4. Shear-driven turbulence does not survive in the heated case, although shooting in the upper direction for C = 4 does still produce a short turbulent transient. In particular, large amplification of the initial disturbance still occurs in the heated case, but the selfsustaining mechanism appears to be disrupted.
To summarise this section, we have observed that a known TW solution of the isothermal pipe flow is suppressed by buoyancy and that it is connected to the transition to turbulence. The observations are consistent with destabilisation of the shear-driven turbulent state, but at this stage another approach is required to forge an approximate quantitative link with the transition from turbulence.

Calculation of the apparent Reynolds number of HHS
In §1.2, where we gave a brief overview of HHS, the (isothermal) equivalent pressure gradient flow (EPG flow) was identified as a useful reference case for heated flows. To calculate the apparent Reynolds number of the EPG reference flow, one must determine the contribution to the mass flux from the buoyancy force that would have been induced in a fixed pressure-gradient flow. Here we summarise the key points of the analysis of HHS and apply it to a selected example case from our data. (The interested reader is referred to sections 3.3 and 3.5 of HHS for a detailed derivation.) In the following section we relate HHS analysis to the phase diagram determined from the simulations of §3.1.
The analysis starts by decomposing the body-force influenced flow (i.e. the total flow) into a pressure-driven flow of equivalent pressure gradient (the EPG reference flow) and  Figure 14: Times series of (left) total dissipation Dtot (normalised by the laminar isothermal value D0 = 2πLz| − 2| = 4πLz) and (right) energy of the streamwise-dependent modes E 3d for simulations started from the lower-branch TW solutions at Re = 3000, α = 1.7 with C = 0 and C = 4. The TW is perturbed by adding ∓ 0.001 (w1 + 0.01w2) (denoted as 'upper' and 'opposite' directions) where w1 and w2 are the first (leading) and second eigenvectors. Shooting in the 'upper' direction leads to turbulence for C = 0, while the flow goes back to laminar when perturbed in the opposite direction. For C = 4 both directions end up at the laminar point. a perturbation flow due to the body force, where the superscripts † and f denote the EPG and the body-force perturbation driven flows, respectively. In contrast to the conventional view, HHS observe that adding a nonuniform (radially-dependent) streamwise body force to a flow initially driven only by a pressure gradient does not alter its turbulent mixing characteristics and in particular the turbulent viscosity remains approximately the same. From this point of view, the bodyforce influenced flow behaves in the same way as the EPG flow and relaminarisation occurs when the Reynolds number Re app of this 'apparent' flow drops below a certain threshold where turbulence cannot be sustained any more. Given the difficulties discussed in §1 to uniquely define a critical Reynolds number for transition, we decided to follow HHS and select a nominal value of 2300, as quoted in many engineering textbooks (see e.g. White 1979). By writing the bulk velocity U b of the EPG flow as the difference between that of the total flow and of the body-force perturbation driven flow, i.e. U † b = 0.5 − U f b , the above relaminarisation criterion can be expressed as (3.2) To determine U f b , the following expression was derived by integrating three times the Reynolds-averaged z-momentum equation of the body-forced perturbation flow: where R f uv (r) := (u z u r ) f is the Reynolds shear stress due to the perturbation flow induced by the body force f (r). The first integral of (3.3), I 1 := 1 2 1 0 (1 − r 2 )f (r) rdr, represents the direct contribution of the body force (which is assisting the flow), while the second integral, I 2 := force perturbed flow is related to that of the total (R uv ) and EPG (R † uv ) flows by using the decomposition (3.1) and is approximated by introducing the eddy viscosity concept, reference flow is calculated using an expression originally suggested by Cess (1958), see (B 2). The resulting eddy viscosity is shown in figure 16(left). By substituting ν † t in (3.6) we can invert for dU f z /dr which plugged into (3.5) gives us the Reynolds stress R f uv (r) (see figure 16(right)). Finally, by inserting the latter in the second integral of (3.3) we obtain Re I 2 = 0.0405. Putting everything together, (3.3) gives U f b = ReI 1 + ReI 2 = 0.12 + 0.0405 ≈ 0.16. Then, using (3.2), Re app = Re 1 − 2 U f b = 2040 < 2300, i.e. the flow is expected to relaminarise. This value obtained for the apparent Reynolds number is reasonable, since relaminarisation occurs after approximately 400 time units (see figure  15 (right)).

HHS prediction of phase diagram and nonlinear dynamics
We now consider the general case of a flow at Re with heating C, while introducing a number of approximations to simplify the analysis.
Firstly, the case discussed in §3.4 (C = 2 and Re = 3000) suggests that Re I 1 has a significantly greater contribution than Re I 2 in determining the body-force perturbation flow. This is found to be generally true for the cases considered herein, as well as those discussed in HHS, and hence we omit the term Re I 2 for simplicity below. The perturbation flow due to the body force can thus be evaluated as where (3.7) has been used for f (r). Secondly, figure 4(bottom right) shows that the temperature mean profiles are remarkably similar in all turbulent shear-driven flows (i.e. ignoring the laminar or convection driven flow states), as far as the integral part of the right-hand side of (3.8) is concerned, despite that the values of the Nu (proportional to the gradient at the wall) are necessarily quite different for different cases. For the case Re = 5300 C = 3.75, for the left-hand side of (3.8) we obtain Re I 1 = 0.164. By applying the above assumption, Let Re app =2300 to find the critical C for flow laminarisation, that is, (3.10) or C cr,1 = 12.5 1 − 2300 Re . (3.11) For C C cr,1 we expect to see rapid transition from the shear-driven turbulent state to the convective state. Noting C = Gr/(16Re), the above can be expressed as a critical Grashof number: Gr cr,1 = 200(Re − 2300) (3.12) Let us now consider the opposite scenario in which the flow under heating C is either laminar or convection driven. Figure 4(bottom right) shows that the temperature profiles in such flows are significantly different from those in a turbulent shear-driven flow, and generally with a much thicker thermal boundary layer, and hence a greater buoyancy force. Consider the extreme case when the radial heat transfer is purely due to conduction and the temperature distribution is given by Θ = r 2 . The buoyancy-driven perturbation flow is therefore Then a second critical C = C cr,2 can be evaluated, below which the flow is expected to transition to the shear-driven turbulent flow. To put it another way, it is predicted that metastability of the shear-driven turbulent state should not be observed for C C cr,2 , so that the turbulent state is stable. Between C cr,1 and C cr,2 the shear-driven state is expected to be metastable, so that this or a convective state may be observed. In terms of the Grashof number, Gr cr,2 = 96(Re − 2300).
(3.15) Equations (3.11) and (3.14) are plotted on the Re − C graph in figure 17 together with all DNS results already presented in figure 7. The data of figure 7 was obtained starting from shear-driven turbulent states. Some additional simulations were performed at Re = 5300 starting from convection-driven states and are reported in figure 17 using hollow symbols, with a slight offset in Re for visualisation reasons. Note that in a Re−Gr graph (3.12) and (3.15) are straight lines (see the inset in figure 17).
Considering a series of DNS runs for a fixed Re, for example Re = 5300, but increasing C values (heating) starting from C = 0, equation (3.11) gives the critical C = C cr,1 above which the flow will be laminarised or switch to convection-driven. On the other hand, starting from a large C when the flow is laminarised or convective, equation (3.14) predicts a critical C = C cr,2 below which the flow will be turbulent when sufficient disturbances are provided in the DNS. As C cr,1 is larger than C cr,2 for a given Re, there is an overlap in the possible state of flow, and consequently there is a hysteresis region in which the flow may or may not be laminarised, depending on the initial flow of the simulation (or experiment). As a result, the Re − C plane can be divided into three regimes by the curves representing the two equations, i.e., turbulent shear-driven flow (regime I), convection-driven or laminar flow (regime III) and regime II in which either of the above may happen dependent on the initial flow. Note that for the Reynolds number range considered here, the linear stability curve (showed as a dashed grey line in figure  7) is always to the right of C cr,2 , i.e. C cr,2 < C LS . The two curves cross at Re ≈ 6000 (not shown), which means that, for Re < 6000 the convective flow is always linearly flow, as in 7, together with equations (3.11) and (3.14) and the linear stability stability curve (dashed red curve in figure 8). Initial conditions are a shear-driven turbulent state, except for the hollow symbols at Re = 5300 which are started with a convection driven state, and similarly cases towards the bottom-right, where it is clear that the shear-driven state decays immediately. stable if C < C cr,2 . Hence, below Re ≈ 6000, shear driven turbulence may be observed for C < C LS .
A plot showing the phase transitions for the fixed Reynolds number Re = 5300 is provided in figure 18, where the Nusselt number is displayed as a function of C for simulations started with either shear-driven or convection-driven states. The two critical C at this Reynolds number, C cr,1 = 7.1 and C cr,2 = 3.4, are indicated with vertical lines in figure 18. Starting from an unheated (C = 0) turbulent flow, applying a low heating (C 7), we observe that the flow remains turbulent over the entire period of simulation (t = 2000). The dynamics thus sits on the upper branch shown in figure  18. As C is increased, the lifetime of shear-turbulence drops below 2000 time units for C 7.5 and turbulence only survives for less than 500 time units at C = 10. It then switches to the convection-type flow. This behaviour is marked in figure 18 by plotting the upper-branch curve with a dashed line for C 7.5 until it crosses the lower-branch at C = 12.5. At this value of C, indeed, the switch to the convective flow appears to be immediate. Now, starting from this convection-driven flow and applying a lower C, the flow remains convection-driven turbulent for C 3.8, or relaminarises for C 3.8. This value of C corresponds to the onset of the linear instability, which is responsible for the kink in Nu as C is decreased. Our previous analysis predicts that for flows on the left of (3.14), their Re app is greater than 2300, hence they may be prone to transition to turbulence subject to sufficient disturbances. Correspondingly, the lower-branch curve in figure 18 is plotted with a dashed line for C < C cr,2 = 3.4 to indicate that in practice (e.g. in a lab experiment) the flow would become shear-driven turbulent again. However, as previously discussed, at this Reynolds number, C cr,2 < C LS . Bistability (between shear or convection driven states) is thus observed for 3.8 C 7.5. The latter value is in very good agreement with the threshold C cr,1 = 7.1 predicted above.  Figure 18: Nusselt number vs C for simulations started with shear and convection initial conditions at Re = 5300. The magenta and cyan vertical lines correspond to the critical buoyancy parameters C cr,1 and C cr,2 given by (3.11) and (3.14), respectively. For values of C C cr1 (C C cr2 ) the shear-driven (convection-driven) state is not supported and correspondingly the upper (lower) branch is plotted with a dashed semi-transparent line.
In figures 19 and 20 the turbulent structures of the isothermal and heated flows at Re = 5300, C = 0 and 5, are compared to those of the EPG reference flow. The latter was computed by performing a DNS with fixed pressure gradient such that Re † p = Re p = 10898.7. The flow structures -streaks and vortices -are visualised as isosurfaces of streamwise velocity and streamwise vorticity fluctuations, normalised by the apparent friction velocity based on the pressure gradient component of the wall shear stress only, u * τ p , where the asterisk * denotes a dimensional quantity here. The resulting apparent friction Reynolds number is Re τ p := u * τ p R * /ν * = Re † τ = 147.6. Comparison between the isothermal and heated flows show that the streaks are relatively unaffected, while vortices are significantly weakened. Our interpretation is that while the streaks are responsible for the saturation of the nonlinearity of the flow, via nonlinear normality of the mean flow (Waleffe 1995), it is relatively 'easy' to produce streaks. Note that the mean axial flow for these cases is almost identical (figure 4), and at the end of §3.3 large initial amplifications of disturbances remains possible in the heated case. It is observed that weaker vortices in the heated case are sufficient to produce saturated streaks of the same amplitude. Thus, vortices are more important in the sense that criticality for transition appears to occur when the vortices are too weak. Comparing now the heated flow with the EPG flow, consistent with the observations of HHS (see their figure 19), it can be seen that the streaks in the heated flow are typically stronger than in the EPG flow, while the vortices are of similar strength. In figure 21 we plot RMS velocity fluctuations. Axial perturbations (a) are not strongly affected by the heating, while the cross-flow components (b) are significantly suppressed. (The plot for u r is very similar to that shown for u θ .) In (c) it is seen that the heated and EPG flow have Figure 19: Three-dimensional visualisations of low (blue) and high (yellow) speed streaks in the isothermal (left), heated (middle) and EPG (right) flows. Isosurfaces of turbulent streamwise velocity normalised by the corresponding apparent friction velocity u z /u τ p = ±4.
very similar cross-components, while axial perturbations in the heated case are slightly stronger than in the EPG flow. These results are consistent with observations from the three-dimensional visualisations of figures 19 and 20, and likewise suggest that it is the weakening of rolls rather than streaks that appear to be responsible for laminarisation.

Conclusions
In this paper we have studied the flow of fluid through a vertically-aligned heated pipe using direct numerical simulations (DNS), linear stability and nonlinear travelling-wave solution analyses. The flow is driven by an externally applied pressure gradient and aided by the buoyancy resulting from the lightening of the fluid close to the heated wall. DNS were performed for a range of Reynolds numbers Re and buoyancy parameters C, where the latter measures the magnitude of the buoyancy force relative to the the pressure gradient of the laminar isothermal shear flow, and three different flow regimes were identified -laminar flow, shear-driven turbulence and convection-driven flow -depending on the flow parameters. At relatively low Re 3500 turbulence is completely suppressed (relaminarised) by buoyancy and as C is increased convection starts driving a relatively quiescent flow. For larger Re, instead, the shear-driven turbulent flow transitions directly to the convection-driven state. Consistent with the appearance of the convective state observed in simulations, a linear instability was found at C ≈ 4, roughly independent of Re for most of the range considered. The result of increasing C can be compared to that of increasing polymer concentration, or Weissenberg number W i, which is known to have a drag reducing effect on turbulent flows (Virk et al. 1967). Similar to our phase diagram (figure 7), a region of relatively quiescent flow has been reported for a certain range of Re and W i (Choueiri et al. 2018;Lopez et al. 2019), although the underlying physical mechanism (elastoinertial instability) is clearly very different from the one studied here (convection driven).
Cases where turbulence is suppressed exhibit a flattened mean streamwise velocity profile. In agreement with recent observations by Kühnen et al. (2018) and Marensi et al. (2019) on the effect of flattening, we found that states that mediate turbulence (lower-branch travelling wave solutions) are "pushed out" from the laminar state, i.e. case correspond closely to its EPG counterpart, while the heated case has slightly stronger streaks.
as C increases, a larger perturbation amplitude or larger Re are required to drive shear turbulence until, for sufficiently large C, the travelling wave is suppressed altogether. Finally, we used the relaminarisation criterion recently proposed by He et al. (2016), based on an "apparent Reynolds number" of the flow, to predict the critical C = C cr,1 (Re) above which the flow will be laminarised or switch to the convection-driven type. This apparent Reynolds number is based on an apparent friction velocity associated with only the pressure force of the flow (i.e. excluding the contribution of the body force/buoyancy). Bistability between shear or convection-driven states was found to occur in the region 4 C C cr,1 where the flow may or may not be laminarised depending on the initial flow of the simulation or experiment.
Comparison of the turbulent flow structures (rolls and streaks) with those of two reference flows -the flow of equivalent pressure gradient (EPG) and that of equivalent mass flux (EFR) -suggests that near criticality for relaminarisation the vortices, rather than the streaks, are more important in the sense that criticality for transition occurs when the vortices are too weak. This picture is not straight forward to reconcile with the interpretation of Kühnen et al. (2018), where relaminarisation is attributed to reduced ability to produce streaks in the presence of the flattened base profile. In the heated case, the base velocity profile does not appear to change significantly while shear-driven turbulece is present. Thus it appears unlikely that transient growth of streaks is affected by the heating. Indeed, laminarisation occurs despite little suppression of the streaks. The experiments of Kühnen et al. (2018) are slightly different, however, in that the var-ious flow manipulations they introduce do change the base profile of the flow. In that case it is correct that transient growth will be affected, although we conjecture that it is the suppression of the vortices due to suppression of the streaks that is responsible for laminarisation in that case. Their numerical experiments in the presence of a force are very similar to the calculations here and of HHS. In that case we expect the mechanism we have described to be more clearly responsible for the laminarisation. gradient. The parameters A + = 27 and κ = 0.42 have been chosen to fit the more recent observations of (McKeon et al. 2005).
For the calculation of §3.4, the (apparent) pressure gradient B and (apparent) Re p are known. The mass flux Re of (B 1) is not yet known, and we wish to determine ν t . An initial estimate for Re is obtained from the approximation of Blasius (1913), which may be written Then, (B 2) can be used to calculate ν t (r), but we must check consistency with (B 1). The latter equation can be inverted for U (r), and, as it has been non-dimensionalised with the same scales of section §2.1, the mean velocity U b = 2 1 0 U (r) r dr should be 0.5. It will not be exactly so, as Re (for the given ∂ z P ) has only been estimated. A better estimate is given by Re := (0.5/U b ) Re, so that ν t can be recalculated and iteratively improved.