2. Numerical method

The solutions are obtained by solving the unsteady three-dimensional incompressible Navier–Stokes equation for a viscous Newtonian fluid:

where $\boldsymbol{v}=(u,v,w)$ and $P$ denote the scaled velocity vector and pressure respectively. Equations (2.1) and (2.2) were non-dimensionalised using the inlet pipe diameter, $d$ , for the length scale and the bulk velocity at the inlet, $U$ , for the velocity scale. The time scale and the pressure scale are therefore $t=d/U$ and ${\it\rho}U^{2}$ , where ${\it\rho}$ is the density of the fluid. The equations are solved with the boundary conditions: corresponding to a fully developed Hagen–Poiseuille flow (2.3) at the inlet, no-slip (2.4) at the walls and an open boundary condition (2.5) at the outlet of the pipe, $\boldsymbol{e}_{z}$ is the unit vector in the axial direction and $r$ is the cylindrical radius. Equation (2.5) enforces Neumann boundary conditions in weak sense for the velocity components, which minimises the possibility of numerical oscillations and reflections of outgoing waves, where $\boldsymbol{n}$ is the normal surface vector directed out of the computational domain. The equations were solved using an open source code nek5000 developed by Fischer et al. (Reference Fischer, Kruse, Mullen, Tufo, Lottes and Kerkemeier2008). Spatial discretisation is based on the spectral-element method using Lagrange polynomials. The equations are reduced to a weak form and discretised in space by Galerkin approximation. Here, $N$ th-order Lagrange polynomial interpolants on Gauss–Lobatto–Legendre points were chosen as the basis for the velocity space, and similarly for the pressure space. In all the simulations $\mathbb{P}_{N}-\mathbb{P}_{N}$ formulations were implemented, which denotes that the same polynomial order was used for both velocity and pressure. The time-stepping in nek5000 is semi-implicit in which the viscous terms of the Navier–Stokes equations are treated implicitly using third-order backward differentiation and the nonlinear terms are treated by a third-order extrapolation scheme (Maday, Patera & Rønquist Reference Maday, Patera and Rønquist1990; Fischer et al. Reference Fischer, Kruse, Mullen, Tufo, Lottes and Kerkemeier2008).

Figure 1. The spectral-element mesh used in the present study with a divergent angle of ${\it\alpha}=26.57^{\circ }$ . (a) Sketch of the domain, (b) cross-section of the mesh (the dark lines represent the elements and the grey lines represent the Gauss–Lobatto–Legendre points) and (c) a three-dimensional view of the mesh near the diverging section. The mesh is made of $K=14\,400$ elements.

Figure 1 shows the geometry of the divergent pipe along with the mesh. It consists of three parts: (i) the inlet, (ii) the diverging section and (iii) the outlet. The velocity field is simulated in the Cartesian coordinate system ( $x,y,z$ ). The expansion ratio is $E=D/d=2$ , where $D$ is the outlet pipe diameter. The length of the divergent section is kept constant in this study and of length $d$ , which leads to a divergence half-angle of ${\it\alpha}=26.57^{\circ }$ . The length of the inlet pipe is $5d$ and the outlet pipe length is $150d$ .

The mesh was developed using hexahedral elements with a non-uniform growth rate. It contains 80 elements with refinement near the wall in the ( $x,y$ ) cross-section and 180 elements in the $z$ direction. A refinement has been applied in the diverging section, as shown in figure 1(c), in order to resolve the flow separation. The streamwise extent of the elements increases along $z$ in the outlet section. The total number of grid points in the simulation is approximately $KN^{3}$ , where $K$ is the number of elements and $N$ is the polynomial order. The flow was initialised with fully developed Poiseuille flow in the inlet section, and each simulation was computed using 512 cores. Table 1 shows the length of the recirculation region for different orders of polynomial at $Re=1000$ . The mesh convergence study was carried out by changing the polynomial order $N$ of the Lagrange polynomial of the spectral elements. The observations used to assess convergence are the flow reattachment point, $z_{r}$ , and the viscous drag, $({\it\rho}U^{2}/2)\mathscr{A}_{w}C_{f}$ , where $\mathscr{A}_{w}$ is the surface area of the outlet pipe wall and $C_{f}$ is the friction coefficient. The length of the recirculation region depends sensitively on the resolution of the separated shear layer, particularly near the separation point. The polynomial order of $N=5$ is sufficient to resolve the flow accurately. These values of $N$ and the mesh have been used in all the following simulations; this corresponds to $KN^{3}\approx 1.8\times 10^{6}$ grid points.

Figure 2. (a) Streamline of the recirculation region of length $z_{r}$ inside the diverging pipe at $Re=300$ , (b) recirculation region length, $z_{r},$ with respect to $Re$ . Here, ● corresponds to best-fit proportionality given by $2z_{r}=0.0866Re$ for the present case, ▫ corresponds to $2z_{r}=0.0874Re$ (Cantwell, Barkley & Blackburn Reference Cantwell, Barkley and Blackburn2010), and ▵ and ▿ correspond to the experimental results (Latornell & Pollard Reference Latornell and Pollard1986; Hammad et al. Reference Hammad, Ötügen and Arik1999)  $2z_{r}=0.088Re$ and $2z_{r}=0.096Re$ respectively for sudden expansions.

To validate the simulations further at higher Reynolds number, the growth of the recirculation region as a function of the Reynolds number is shown in figure 2. The simulations show that the extent of the recirculation region is of the form $2z_{r}=0.0866Re$ , which agrees well with previous studies for sudden-expansion flow. Unlike the sudden-expansion flow, for a divergent pipe, the recirculation is formed after a critical flow velocity, shown in the inset to figure 2, which depends upon the divergence half-angle ${\it\alpha}$ and $Re$ (Peixinho & Besnard Reference Peixinho and Besnard2013).

3. Results and discussion

For sudden-expansion pipe flow, numerical simulations (Cantwell et al. Reference Cantwell, Barkley and Blackburn2010; Sanmiguel-Rojas et al. Reference Sanmiguel-Rojas, Del Pino and Gutiérrez-Montes2010) have shown that the flow is unstable to infinitesimal perturbations for $Re\approx 3273$ , but the transition in experiments occurs at much lower $Re$ (Sreenivasan & Strykowski Reference Sreenivasan and Strykowski1983; Latornell & Pollard Reference Latornell and Pollard1986; Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009). The exact nature of the observed instability is therefore unclear. Small disturbances in an experimental set-up are likely to be amplified due to the convective instability mechanism, and appear to be necessary to realise time-dependent solutions. Numerical results (Cantwell et al. Reference Cantwell, Barkley and Blackburn2010) have shown that small perturbations are amplified by transient growth in the sudden expansion for $Re\leqslant 1200$ , advect downstream and decay. Here, the initial simulations showed that the flow is linearly stable for up to $Re\gtrsim 2200$ for the present computational domain. For larger $Re$ , the recirculation bubble extends close to the end of the outlet section and cannot be calculated reliably. In order to induce early transition, a disturbance is applied to the numerical system in the form (Sanmiguel-Rojas & Mullin Reference Sanmiguel-Rojas and Mullin2012)

adding a finite-amplitude crosswise velocity of magnitude ${\it\delta}$ .

The perturbation (3.1) distorts the flow, nudging it towards the $y$ direction. A perturbation value of ${\it\delta}=0.001$ is applied in the following simulations. For the sudden-expansion pipe, this is the value of ${\it\delta}$ for which most results are presented in Sanmiguel-Rojas & Mullin (Reference Sanmiguel-Rojas and Mullin2012). Results were found to be compatible with the imperfections found in experiments. Figure 3(a–c) shows cross-sections of the pipe at $z=22.5$ and presents contours of the axial velocity. Figure 3(a) is at $Re=1000$ where the flow remains almost axisymmetric. For $Re=1600$ , figure 3(b) shows an asymmetry which can just be identified by comparing the solid and dashed lines. To see more clearly the effect of perturbation on the flow, the perturbed flow is subtracted from the (unperturbed) base flow, where it can be observed in the contour plot, figure 3(c), that the flow is accelerated on the right-hand side of the pipe section and decelerated on the opposite side. The applied perturbation at the inlet creates a recirculation region with a biased extent (see figure 3 d). The reattachment pattern is very sensitive to the form of the perturbation given at the inlet, which motivates the application of a simple form of disturbance.

Figure 3. Cross-sections of the pipe. Contour line plots of the axial velocity (solid black lines) taken at $z=22.5$ for (a) $Re=1000$ and (b) $Re=1600$ . The dashed line corresponds to the inlet pipe diameter and the blue lines with arrows represent the crosswise velocities within the recirculation region. (c) Contour plot of the perturbation, i.e. flow with perturbation ( ${\it\delta}=0.001$ ) subtracted from the base flow ( ${\it\delta}=0$ ) for $Re=1600$ . (d) Streamwise cross-section of the flow around the reattachment point at $Re=1600$ with ${\it\delta}=0.001$ .

Figure 4. (a) Asymmetry growth of the flow measured by the square of the distance of the centroid from the centre of the pipe, ${\it\varepsilon}$ , as a function of $Re$ . The lines are least-squares fits of the data and the intersection of the lines is at $Re_{c}=912$ for the estimate of the symmetry-breaking bifurcation point. (b) Oscillations of ${\it\varepsilon}$ at $Re=1650$ as a function of time. The inset is the fast Fourier transform of the signal with a fundamental frequency $f=0.468$ and a period doubling sub-harmonic $SH=0.234$ .

The asymmetry growth of the flow in the cross-section at $z=22.5$ is measured by calculating the distance of the position of the peak axial velocity component from the centre of the pipe. The square of this distance is denoted ${\it\varepsilon}$ (Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009). Figure 4(a) shows ${\it\varepsilon}$ as a function of $Re$ with least-squares fitting on the data obtained. It can be seen that at low $Re$ there is no variation in the position of the centroid; a steady symmetric state is observed for $Re<912$ . As $Re$ increases, a symmetry-breaking bifurcation occurs at a critical $Re_{c}=912$ . This value is smaller than the case of a sudden-expansion pipe ( $Re_{c}=1139$ in the experiment by Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009). Clearly, the critical $Re$ depends on ${\it\alpha}$ and ${\it\delta}$ . The value of ${\it\varepsilon}$ increases linearly $(912<Re<1500)$ , forming a steady asymmetric state, with biased growth in the recirculation region. The magnitude of the symmetry deviation grows as the square root of $Re$ , typical of supercritical bifurcation. At larger $Re$ an oscillation state arises ( $1500\leqslant Re\leqslant 1650$ ), and the flow becomes time-dependent, due to the spatio-temporal oscillation of the reattachment point downstream (Sreenivasan & Strykowski Reference Sreenivasan and Strykowski1983). The error bars in figure 4(a) represent the amplitude of the fluctuations in ${\it\varepsilon}$ . These oscillations are also observed in experiments on sudden-expansion flow (Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009). When $Re=1650$ , the flow experiences quasi-periodic oscillations in the shear layer around the recirculation region. This can be seen in the velocity components along the axial as well as in the crosswise direction, and also in the ${\it\varepsilon}$ evolution (see figure 4 b). A fast Fourier transform (FFT) was performed on the signal to identify the dominant frequency. The inset in figure 4(b) is the FFT of the signal as a function of frequency, where $f=0.468$ and $SH=0.234$ is a period doubling sub-harmonic. Here, $f$ seems to correspond to the frequency of vortex shedding around a circular or spherical body which occurs due to the Kelvin–Helmholtz instability (Fabre, Auguste & Magnaudet Reference Fabre, Auguste and Magnaudet2008; Bobinski, Goujon-Durand & Wesfreid Reference Bobinski, Goujon-Durand and Wesfreid2014). This frequency of oscillation depends upon the type of the perturbation added to the system (Marquet et al. Reference Marquet, Sipp, Chomaz and Jacquin2008; Ehrenstein & Gallaire Reference Ehrenstein and Gallaire2009).

Figure 5. (a) Plot of the friction coefficient, $C_{f}$ , with respect to $Re$ . The blue filled circles represent the steady laminar asymmetric flow and the red filled squares represent the unsteady localised turbulent state. The continuous line represents the fit for the laminar state, $C_{f}=1.97/Re-0.0012$ , and the dotted line represents a fit for the localised turbulent state, $C_{f}=0.0066(Re\times E)^{-0.22}$ . The shaded region $1650\lesssim Re\lesssim 3273$ is the coexistence regime (hysteresis for ${\it\delta}=0$ ), where the left-hand subregion indicates the extent of the regime explored on the laminar branch in the present system. (b) Contour plot of the streamwise velocity of the localised turbulence at $Re=1680$ with ${\it\delta}=0.001$ .

Figure 5(a) shows the friction coefficient, $C_{f}$ , as a function of $Re$ , computed on the wall of the outlet section. At low $Re$ , the flow is steady and asymmetric, and the value of $C_{f}$ decreases. A significant contributor to the low values of $C_{f}$ is that the present flow includes the recirculation region, which extends up to approximately half of the outlet section before transition. Around the transition Reynolds number, $Re_{t}\simeq 1680$ , the recirculation region inside the pipe breaks and leads to a localised turbulent state, shown in figure 5(b).

An important feature of the turbulence observed here is that it remains spatially localised at a constant position, as observed in sudden-expansion pipe flows (Sreenivasan & Strykowski Reference Sreenivasan and Strykowski1983; Sanmiguel-Rojas & Mullin Reference Sanmiguel-Rojas and Mullin2012). The formation of turbulence, near the diverging section, increases $C_{f}$ due to higher internal mixing and resulting shear at the boundary. In this regime, $C_{f}$ values scale roughly with the same exponent as the Blasius friction law, even though straight pipe flows are not turbulent at these $Re$ values. The present calculations for the perturbed flow were run up to $t=600$ and the localised turbulence remained present. The perturbation was then removed and the flow was simulated up to $t=1200$ . The turbulence was observed to be self-sustained, and to occupy the same spatial position. The left-hand shaded region in figure 5(a) shows the range of $Re$ in which a laminar state as well as a turbulent state coexist for $1650\lesssim Re\lesssim 2200$ for the present computational domain. Simulation above $Re>2200$ without perturbation produces a steady laminar flow with a recirculation region that extends close to or beyond the outlet section. We have therefore limited the range of $Re$ for computation on the laminar branch. For the case of simulations with a perturbation, the amplified energy in the diverging section breaks the recirculation region, creating an early transition, forming localised turbulence, and the computation may be carried out for larger $Re$ along the turbulent branch. Global stability analysis (Sanmiguel-Rojas et al. Reference Sanmiguel-Rojas, Del Pino and Gutiérrez-Montes2010) has revealed that the first bifurcation for the sudden-expansion pipe occurs at $Re\gtrsim 3273$ , above which natural transition can be expected without any added perturbation. Given the much larger computational cost, and that we have already computed a range of $Re$ where the laminar and turbulent flows coexist, we have not pursued the linear instability.

Figure 6. Relaminarisation study. (a) Spatio-temporal diagram of the streamwise vorticity along the centreline of the pipe, where $z=0$ corresponds to the start of the diverging section. The localised turbulence decays from $Re_{0}=2000$ to $Re_{f}=1500$ . (b) Relaminarisation time, $t_{R}$ , versus $Re_{f}$ .

Figure 7. The evolution of localised turbulence. Spatio-temporal diagram of the streamwise vorticity along the centreline of the pipe, where $z=0$ corresponds to the start of the diverging section, for (a) $Re=4000$ , (b) $Re=3000$ , (c) $Re=2000$ and (d) $Re=5000$ . Panels ( $e$ ) and (  $f$ ) show $q=\sqrt{v^{2}+w^{2}}$ in (red) and the streamwise velocity $u$ (blue) at the final time step of (c) and (d). The dashed lines represent the Poiseuille centreline velocity in the inlet and outlet sections.

Further relaminarisation simulations were performed, where localised turbulence was generated at $Re_{0}=2000$ and the decay to laminar flow was observed for $Re$ below $Re_{t}$ . In figure 6(a), the spatio-temporal diagram shows a typical relaminarisation case. At $t=0$ , $Re$ is reduced suddenly from $Re_{0}=2000$ to $Re_{f}=1500$ . Here, the localised turbulence detaches from the inlet section almost immediately, then convects downstream and simultaneously decays, which can be seen as the disappearance of the vortical structures (Sibulkin Reference Sibulkin1962; Sreenivasan Reference Sreenivasan1982). The relaminarisation time, $t_{R}$ , was obtained by monitoring the time taken for the total energy in the computational domain to fall below a threshold of $10^{-6}$ . Above $Re_{f}=1500$ the turbulence leaves the computational domain before falling below the threshold. Figure 6(b) shows $t_{R}$ as a function of $Re_{f}$ . The straight line fit indicates that the decay time of the turbulence increases linearly for $Re<1500$ , as identified in experiments (Peixinho & Besnard Reference Peixinho and Besnard2013). Here, no significant period of time was observed before the detachment of turbulence from the walls. Simulations were not carried out within the hysteresis region due to high computational cost. For these $Re$ values the turbulence is self-sustained for some time before detachment from the inlet section, and the time before detachment is expected to diverge rapidly as in uniform pipe flow (Avila et al. Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011).

Figure 7 shows the streamwise vorticity spatio-temporal diagram of simulations for (a) $Re=4000$ , (b) $Re=3000$ , (c) $Re=2000$ and (d) $Re=5000$ ; the horizontal axis represents the centre axis of the pipe from the diverging section to the outlet. The streamwise vorticity value has been normalised with the maximum vorticity and been plotted with the same scale for comparison purposes. It can be seen that for $Re=2000$ the turbulence onsets at $t=25$ and initially moves downstream; at $t=100$ the turbulence starts moving upstream towards the diverging section and finally holds a stable position $z\simeq 10$ . For $Re=5000$ , the onset of turbulence occurs at nearly the same time as that of $Re=2000$ , but the amount of time it takes to reach a localised position is $t=40$ , which is much smaller than that for $Re=2000$ . The time taken to hold a stable position decreases as $Re$ increases. The velocity trace downstream of the localised turbulence for $Re=2000$ recovers laminar flow (see figure 7 e). It should be noted that the Reynolds number based on the outlet diameter is half the value of $Re$ . In the case of $Re=5000$ , the flow downstream of the intense region of turbulence exhibits small patches of intense vorticity (see figure 7 d). The streamwise velocity trace (see figure 7 f) suggests weak turbulence, which does not return to laminar flow and eventually could lead to puff splitting (Avila et al. Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011; Shimizu et al. Reference Shimizu, Manneville, Duguet and Kawahara2014). This property of expansion flow with laminar inlet profile forming localised turbulence and decaying in the outlet section is in good agreement with experiments (Peixinho & Besnard Reference Peixinho and Besnard2013).

Figure 8. Spacial correlation on the streamwise velocity of localised turbulence at $Re=2000$ . Contour of correlation, $C({\it\theta},z)$ , at (a) $r=0.5d$ and (c) $r=0.8d$ . Projection function of $C({\it\theta},z)$ for different azimuthal wavenumbers $m$ at (b) $r=0.5d$ and (d) $r=0.8d$ . (e–h) Cross-sections of the axial flow relative to the time averaged profile with fast flow (white) and slow flow (red) taken at the corresponding vertical dashed lines.

The structure within the localised turbulence is further studied using spatial correlations, which have been used to identify fast and slow streaks that dominate the coherent structures within puffs in pipe flow (Willis & Kerswell Reference Willis and Kerswell2008). The correlation in the streamwise velocity is obtained using the function

where $\langle \cdot \rangle _{s}$ indicates averaging over the subscripted variable, $u_{z}$ is the instantaneous axial flow velocity and $r$ is the radial position. The signature of structures of a particular azimuthal wavenumber $m$ is obtained by projecting the correlation function, $C_{m}(z)=2\langle C({\it\theta},z)\cos (m{\it\theta})\rangle _{{\it\theta}}$ (Willis & Kerswell Reference Willis and Kerswell2008). Figure 8(a,b) shows the correlation at $r=0.5d$ , and it can be seen that the $m=1$ mode dominates the flow, whereas in figure 8(c,d) at $r=0.8d$ , the $m=2$ structure dominates the flow along with $m=3$ with a much smaller correlation value, which suggests that the flow is more active in the centre region than near the wall. Overall this analysis points out that the localised turbulence in the gradual expansion possesses a different flow structure from the turbulent puff (Wygnanski & Champagne Reference Wygnanski and Champagne1973; Willis & Kerswell Reference Willis and Kerswell2008; Shimizu & Kida Reference Shimizu and Kida2009), where $m=3$ and 4 dominate the flow near the wall. The cross-sections in figure 8(e–h) indicate slow and fast moving flow.