1. Introduction
Pipes with an expansion (figure 1) have flow features shared by many engineering and biomedical applications, with recirculation and reattachment regions, and transition between laminar and turbulent flow. Typical engineering applications cited include ramjets and combustors (Sreenivasan & Strykowski Reference Sreenivasan and Strykowski1983) and ventilation systems (Lukács & Vad Reference Lukács and Vad2021). Related biomedical examples include blood flow in arteries (Zmijanovic et al. Reference Zmijanovic, Mendez, Moureau and Nicoud2017) and air flow in the trachea (Brouns et al. Reference Brouns, Jayaraju, Lacor, De Mey, Noppen, Vincken and Verbanck2007).
Axisymmetric schematic of the flow in a pipe with a sudden expansion.
The turbulence puffs observed in pipes with a sudden expansion are different from those in transitional flow in straight pipes, where the puff is formed then moves downstream with the flow (Nishi et al. Reference Nishi, Ünsal, Durst and Biswas2008), whereas in the pipe with an expansion, the puff can stabilize and persist at a certain location downstream of the expansion, while still being surrounded by laminar flow upstream and downstream.
Historically, experimental studies in pipes with an expansion (Back & Roschke Reference Back and Roschke1972; Sreenivasan & Strykowski Reference Sreenivasan and Strykowski1983; Latornell & Pollard Reference Latornell and Pollard1986; Pak, Cho & Choi Reference Pak, Cho and Choi1990; Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009; Lebon et al. Reference Lebon, Peixinho, Ishizaka and Tasaka2018b) have shown that, as flow rates are increased to inlet Reynolds numbers 
$Re$ of the order of 1000, the flow becomes asymmetric then breaks down into an unsteady state, which is characterized by transition from a steady symmetric laminar state at the inlet, to an unsteady asymmetric turbulent state some distance after the expansion.
On the other hand, numerical studies (Sanmiguel-Rojas, del Pino & Gutiérrez-Montes Reference Sanmiguel-Rojas, del Pino and Gutiérrez-Montes2010; Cliffe et al. Reference Cliffe, Hall, Houston, Phipps and Salinger2012; Sanmiguel-Rojas & Mullin Reference Sanmiguel-Rojas and Mullin2012; Selvam, Peixinho & Willis Reference Selvam, Peixinho and Willis2015, Reference Selvam, Peixinho and Willis2016; Nguyen et al. Reference Nguyen, Shadloo, Hadjadj, Lebon and Peixinho2019) have typically required relatively higher flow rates, with 
$Re$ closer to 2000 and up to 5000, and/or the addition of numerical inlet perturbations in order to observe the transition from laminar to turbulent flow. In addition, the position where the asymmetry begins also occurs overall further downstream compared to the physical experiments. The earlier occurrence of turbulence in the physical experiments, in terms of both flow rate as well as position in the pipe downstream of the expansion, can be explained by the fact that the numerical simulations are free of any so-called imperfections or disturbances, which inevitably exist in the real world. These can be small pipe wall scratches or defects, flow disturbances caused by the pump that moves the fluid, or even vibrations from (un)related adjacent equipment. In numerical simulations, the main disturbances observed are arbitrarily added perturbations (other than those arising from numerical calculation errors), which are expected to represent to some degree these real-world imperfections.
In the straight-pipe literature (i.e. without an expansion), the seminal experimental study of Reynolds (Reference Reynolds1883) investigated the transition between the laminar and turbulent states, and observed the remarkable sensitivity of transition to disturbances. More recently, the studies of Meseguer & Trefethen (Reference Meseguer and Trefethen2003), Hof et al. (Reference Hof, Westerweel, Schneider and Eckhardt2006) and Peixinho & Mullin (Reference Peixinho and Mullin2007), as well as the review of Eckhardt et al. (Reference Eckhardt, Schneider, Hof and Westerweel2007), lead to the notion that, unless the numerical simulation of a pipe is perturbed ‘enough’, it will not develop from the laminar state into a turbulent one, regardless of the flow rate.
Conversely, the numerical simulations of Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018), for example, did not require the addition of numerical inlet perturbations in order to observe turbulence. Other numerical studies that did employ perturbations also found a critical or natural transition point where, regardless of the manual addition of perturbations, increasing the flow rate caused turbulence (Sanmiguel-Rojas et al. Reference Sanmiguel-Rojas, del Pino and Gutiérrez-Montes2010; Sanmiguel-Rojas & Mullin Reference Sanmiguel-Rojas and Mullin2012; Selvam et al. Reference Selvam, Peixinho and Willis2015, Reference Selvam, Peixinho and Willis2016). To the authors’ knowledge, this behaviour has not been explained, except for the observation that there are background disturbances and numerical noise that could cause this. Possible sources of disturbance other than the inlet boundary condition could have been numerical errors, such as truncation and aliasing errors, arising from the computational grid/mesh or from the numerical procedures themselves. Such errors are expected always to be present in numerical simulations of turbulence. More discussion about truncation and aliasing errors can be found in Ghosal (Reference Ghosal1996) and Kravchenko & Moin (Reference Kravchenko and Moin1997). This leads to the question of what, and to what extent, different numerical parameters and methods in the simulation influence the onset of turbulence. In this context, the present paper reports the results of an investigation of the flow in a pipe with an expansion that specifically explores sources of numerical perturbations in the simulations that can lead to transition to turbulence.
In this study, the numerical settings that were found to significantly influence the transition across different flow rates were the time step and the grid sizes, which are related, respectively, to the temporal and spatial resolution of computational fluid dynamics (CFD) simulations. The influence of these as sources of numerical error and perturbation are evaluated and compared against simple numerical perturbations added to the inlet boundary condition. Centreline streamwise velocity and velocity fluctuations are evaluated, in addition to reattachment lengths and the turbulence kinetic energy budget of the puff characteristic of this flow. The results are also compared to literature experimental and numerical studies.
2. Numerical methods and modelling
The simulations were run with a laminar Poiseuille flow inlet at bulk velocity (
$U$) equivalent to Reynolds numbers varying from 100 to 2500 (
$Re = U d / \nu$) with a pipe axisymmetric expansion ratio 
$D/d = 2$ (see figure 1). The streamwise lengths of the computational domain used were 
$L_i = 10 h$ and 
$L_o = 200 h$, with further discussion on these values presented in § 4. A three-dimensional (3-D) representation of the geometry with the centreline along which variables were evaluated can be seen in figure 2.
Grid G3: (a) frontal view of the inlet region shaded in darker colour and (b) side view of the region around the sudden expansion. (c) The 3-D geometry and location of the centreline and the origin of the coordinate axes.
The simulations were run using the CFD toolbox ‘OpenFOAM’, which is the most widely used open-source CFD solver available, and has already been used for direct numerical simulations (DNS) of a pipe with a sudden expansion, for example by Sanmiguel-Rojas & Mullin (Reference Sanmiguel-Rojas and Mullin2012) and Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018). This solver is based on the flux-conservative discretization method known as the finite-volume method, which is used to approximate the Navier–Stokes partial differential equations as a matrix of algebraic equations, through the divergence or Gauss’ theorem. Using Gibbs’ dyadic notation, the incompressible continuity and Navier–Stokes equations are, respectively,
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \boldsymbol{u} \boldsymbol{u} ) - \boldsymbol{\nabla} \boldsymbol{\cdot} ( \nu \boldsymbol{\nabla} \boldsymbol{u} ) ={-} \boldsymbol{\nabla} p , \end{gather}
where 
$t$ represents the time, 
$\boldsymbol {u}$ is the instantaneous velocity vector, 
$\nu$ is the kinematic viscosity and 
$p$ is the kinematic pressure. No turbulence closure modelling was added, but instead the DNS approach was used for the present study.
The fluid domain was discretized into a computational grid composed of hexahedral cells using ANSYS Meshing v194, with frontal and side views in figure 2. In total, six different grids were used in the present study, with up to 14.6 million control volumes or cells. Grid details are available in table 1, with the cylindrical coordinate system (
$r,\theta,x$) transformed from the Cartesian system (
$x,y,z$) using 
$x$ as the cylindrical axis, 
$r=\sqrt {y^2+z^2}$ and 
$\theta = \arctan (y/z)$. The distribution of the grid cell sizes in the radial 
$r$ and streamwise 
$x$ directions can be seen in figure 3. When changing grid refinement, the grids were refined or coarsened uniformly across the entire domain, maintaining the same aspect ratios and element-to-element size ratios shown in figure 3. Aspect ratios were below 10 and maintained across all grids, with ratios within 1 to 4 starting from the sudden expansion region increasing up to 3 to 10 at the outlet boundary. The highest aspect ratios correspond to where the elements are thinner in the 
$r$ direction, 
$r \rightarrow R/2$ and 
$r \rightarrow R$, as seen in figure 3(a).
Distribution of grid element sizes for grid G1 in the (a) radial and (b) streamwise directions; 
$\mathrm {d}r$ and 
$\mathrm {d}\kern0.06em x$ indicate the element sizes in each respective direction, and 
$R=D/2$.
Total number of control volumes and number of control volumes in the cylindrical 
$r$, 
$\theta$ and 
$x$ directions, before and after the sudden step expansion, for each grid.
Over a hundred simulations were run for this study, and the simulations with the finest grids took 2 to 3 months of computational time each using one compute node with two Intel Xeon Gold 6148 processors, a total of 40 CPU physical cores in parallel in a single node, available on the University of Saskatchewan's Plato high-performance computing cluster. The version 7 of OpenFOAM was compiled using GCC 5.4.0, with solver parallelization through OpenMPI 2.1.1.
The PISO algorithm, originally proposed by Issa (Reference Issa1986), was used to solve the system of equations. This is a non-iterative algorithm, adjusting pressure–velocity as a whole at each stage to satisfy continuity and momentum simultaneously. The algorithm consists of performing one predictor stage, followed by two or more corrector stages for each time step, without additional outer iterations, without relaxation factors and without gradient limiters, so that the solution of each time step is a legitimate approximation to the differential equations. According to Issa (Reference Issa1986), using more than two corrector stages is expected to result in higher than second-order temporal accuracy for the solution of the pressure and velocity fields. After preliminary tests at 
$Re = 2500$ evaluating the computational costs and precision of the simulations, a good compromise favouring precision was found with five corrector stages (compared to two) and scaled residual convergence criteria of 
$1\times 10^{-6}$ (compared to 
$1\times 10^{-8}$) for pressure and velocity.
Implicit second-order backward differencing was used to approximate the temporal derivatives and second-order central differencing was used for the spatial derivatives. The OpenFOAM generalized geometric–algebraic V-cycle multigrid linear solver with a Gauss–Seidel smoother was used to solve the system of equations.
In order to avoid negatively influencing the consistency of the results, the initial condition for all simulations at 
$t=0$ for the fluid domain was the same, with the fluid at rest with zero kinematic pressure and zero velocity. No-slip boundary conditions were applied at the walls. The inlet boundary condition was the steady-state Poiseuille velocity profile (with its characteristic shape schematized in figure 1), which is the theoretical solution for steady-state laminar flow in a straight pipe, with flow rates varying from Reynolds numbers 100 to 2500, as mentioned previously. At the outlet surface, an advective or convective outflow boundary condition is applied, which is based on solving the equation
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + ( \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{u} = 0 , \end{equation}
where 
$\boldsymbol {u}_{n}$ is the advection velocity calculated from the flow rate normal to the outlet boundary. This type of boundary condition, similar to the one used by Vittal Shenoy et al. (Reference Vittal Shenoy, Safdari Shadloo, Peixinho and Hadjadj2019), with further discussion in Ol'shanskii & Staroverov (Reference Ol'shanskii and Staroverov2000) and Hasan, Anwer & Sanghi (Reference Hasan, Anwer and Sanghi2005), is used to minimize numerical oscillations or reflections of outgoing waves and vortices, but is otherwise expected to behave the same as a zero-velocity-gradient condition. The kinematic pressure at the outlet is calculated from the velocity considering the dynamic pressure (i.e. 
$| \boldsymbol {u} | ^2 / 2$).
When indicated, pseudo-random fluctuations with limited amplitudes were superimposed onto the inlet Poiseuille velocity profile as a simple perturbation in every direction, updated every time step, as follows:
\begin{gather} \boldsymbol{u} = \boldsymbol{u}_{inlet} + \boldsymbol{\delta} | \boldsymbol{u}_{inlet} | , \end{gather}
\begin{gather}\delta_i = r_i \mathcal{A} . \end{gather}
Here 
$\boldsymbol {u}_{inlet}$ is the specified inlet Poiseuille velocity vector, 
$\boldsymbol {\delta }$ is the perturbation vector and 
$r_i$ indicates randomly generated numbers uniformly distributed in the interval from 
$-1$ to 
$+1$. The idea behind these fluctuations is for them to represent an effect comparable to turbulence at the inlet, in a procedure similar to the one described by Zmijanovic et al. (Reference Zmijanovic, Mendez, Moureau and Nicoud2017). Knowing that the root-mean-square (r.m.s.) of the velocity fluctuations caused by this distribution is equal to the standard deviation of the distribution, an amplitude of 
$\mathcal {A}=1\,\%$ (percentage of the mean velocity, from 
$-1\,\%$ to 
$+1\,\%$) will result in a turbulence intensity of 
${\approx }0.6\,\%$ (r.m.s. velocity fluctuations divided by the mean velocity). Throughout the simulations, perturbation amplitudes 
$0.001\,\% \leqslant \mathcal {A} \leqslant 50\,\%$ were used for different 
$Re$. An example visualization of these perturbations in the inlet velocity profile is given in figure 4. Note that these perturbations can change the flow rate instantaneously, but will not change the average flow rate.
Inlet boundary condition with Poiseuille profile and example shaded region indicating range of possible fluctuations; instantaneous results sampled from the simulation with 
$Re = 1500$ and 
$\mathcal {A}=5\,\%$ on grid G4, and 
$R_i = d/2$.
The time-step sizes (
$\Delta t$) are the same for the entire grid and represented using the Courant–Friedrichs–Lewy (CFL) number to non-dimensionalize them, which is calculated, for each control volume or cell in the grid, as
\begin{equation} \textrm{CFL} = \frac{\Delta t \sum_f |\boldsymbol{u}_f \boldsymbol{\cdot} \boldsymbol{S}_f|}{2 V} , \end{equation}
where 
$\sum _f$ indicates the sum over all faces, 
$\boldsymbol {u}_f$ is the face velocity vector, 
$\boldsymbol {S}_f$ is the face surface vector area and 
$V$ is the cell volume. The maximum CFL value across the grid for each simulation is reported here.
Similar to Wu & Moin (Reference Wu and Moin2008) and Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018), the start-up impact of the unrealistic initial condition was minimized by starting the simulations with time steps smaller than maximum CFL numbers of 0.05, and then slowly increasing them in changes smaller than 20 % until hitting the simulations’ targeted maximum CFL.
The flow-through time 
$\tau = U_o/L_o$ was calculated based on the bulk velocity 
$U_o$ along the length of the pipe after the expansion 
$L_o$; see figure 1. The simulations were run for at least one flow-through time before taking time averages, and the time averages used the data from at least another flow-through simulation time. Time-averaged results in this study therefore are from at least 
$t/\tau = 2$, using data starting from 
$t/\tau = 1$. Instantaneous results were sampled at 
$t/\tau =2$, unless otherwise specified. For reference, in this geometry, one flow-through time based on 
$L_o$ is equivalent to 50 flow-through times based on the outlet pipe diameter 
$D$. These times were found to be sufficient to observe stable positioning of the turbulence puff or the development of fully laminar flow (see e.g. figure 8 later). Some tests were run with time averaging up to 
$t/\tau = 10$, but no significant changes compared to 
$t/\tau = 2$ were observed.
3. Results and discussion
Two main numerical settings were found to be capable of causing turbulent behaviour other than the addition of numerical perturbations at the inlet: grid and time-step sizes. These are respectively associated with the spatial and temporal discretization procedures that are inherent to time-resolved CFD simulations. To the authors’ knowledge, the influence of these parameters as perturbations on the turbulent behaviour of DNS has not yet been documented. This section starts by discussing the general flow features observed in this flow, followed by the impact of grid and time-step sizes on the turbulent behaviour of the simulations, then examining the turbulence kinetic energy budget predicted by the simulations. Finally, the effects of numerical inlet perturbations as well as a comparison to literature results is presented.
3.1. Flow features
At the tested Reynolds numbers of 100 to 2500, without additional perturbations the simulations would naturally result in the laminar steady solution as seen in figure 5(a), whereas with sufficient perturbations they would result in a turbulent unsteady solution as seen in figure 5(b). In the turbulent solution, the flow transitions from laminar flow to turbulence at a position downstream from the sudden expansion, with the formation of a turbulence cloud or puff. After the puff, the flow transitions back to laminar flow again. This turbulence puff is stable at a position downstream of the expansion that does not change significantly with time. Using the definitions for straight-pipe flow discussed by Barkley (Reference Barkley2016) and Song et al. (Reference Song, Barkley, Hof and Avila2017), this puff has a strong upstream front and unstable dynamics in the sense that the flow has turbulent and laminar states that are intermittent in position (i.e. there is the transition from laminar flow to turbulence then back to laminar as the flow moves downstream). The flow with the turbulence puff in the pipe with an expansion has a similar behaviour to the turbulence puff and to the turbulence slug at low 
$Re$ in straight pipes, as these have a weak downstream front, with subsequent relaminarization of the flow. Still in comparison to the turbulence puff in a straight pipe, the main difference is that the puff in the pipe with an expansion is localized and does not move with the flow, whereas the puff in the straight pipe moves downstream with the flow.
Instantaneous velocity contour results (a) without and (b) with inlet perturbations of 
$\mathcal {A} = 0.7\,\%$, using grid G4 and 
${\rm CFL} = 0.5$ at 
$Re= 2500$.
In figure 6 it is possible to visualize, through the 
$\lambda _2$ criterion (Jeong & Hussain Reference Jeong and Hussain1995), the core of the vortices produced by this flow, coloured by the vorticity magnitude. There is the formation of a toroidal or ring-like vortex immediately downstream of the sudden expansion, which is observed regardless of turbulence or the presence of perturbations at the inlet. This expansion vortex can also be observed in the results of Vittal Shenoy et al. (Reference Vittal Shenoy, Safdari Shadloo, Peixinho and Hadjadj2019). With inlet perturbations, in the laminar region right after the expansion there is also the formation of a larger clockwise spiral or helical vortex that grows until the flow breaks down into the turbulence puff with numerous smaller vortices. Further downstream from the puff, the vortices elongate and assume shapes similar to hairpin or horseshoe vortices, while also reducing in number as the flow relaminarizes.
Visualization of the cores of vortices at 
$\lambda _2 = -0.5\ {\rm s}^{-2}$ coloured by the vorticity magnitude (a) without and (b) with inlet perturbations of 
$\mathcal {A}= 0.7\,\%$ using the finest grid G1 and maximum 
${\rm CFL} = 0.5$ at 
$Re =2500$.
Figure 7 shows more clearly the region around which the recirculation region ends (i.e. where the flow reattaches), through velocity contour lines. At 
$Re=800$, increased perturbations and subsequent transition to turbulence did not change the streamwise length of the recirculation region significantly, as seen in figure 7(a) for the laminar solution compared to figure 7(b), where the inlet perturbations were sufficient to trigger turbulence. Whereas in the laminar flow the recirculation region is thin and slowly tapers off, in the turbulent flow the recirculation region is thicker and ends abruptly. Additional discussion on figure 7, in the context of other simulations, is presented in § 3.5. Lastly, figure 7(b) also depicts contours of the same large laminar clockwise helical vortex in the recirculation region that was observed in figure 6(b).
Instantaneous streamwise velocity 
$u_x$ contour lines with inlet perturbation amplitudes of (a) 
$\mathcal {A}= 5\,\%$ and (b) 
$\mathcal {A}= 20\,\%$, using grid G4 and maximum 
${\rm CFL} = 0.5$ at 
$Re = 800$. Blue lines correspond to velocity below, red above and black equal to zero. This sectional view is from around 
$x/h \approx 45$ to 
$75$, with the 
$y$ dimension stretched to 150 %.
Figure 8 shows instantaneous snapshots of centreline streamwise velocity plotted as a function of the position 
$x/h$ downstream of the expansion, with the time 
$t$ non-dimensionalized by the flow-through time 
$\tau$, for different inlet perturbation amplitudes. Figure 8(a) shows the evolution of a simulation with a perturbation amplitude 
$\mathcal {A} = 0.01\,\%$ that is not enough to sustain turbulence: initially, a turbulence puff is formed near the expansion due to the start-up of this flow (as seen by the sharp drop as well as the oscillations in the centreline velocity), but the puff is advected downstream as the flow laminarizes behind it. After approximately one flow-through time, the flow is fully laminar, resulting in a reattachment length of 
$L_R/h \approx 133$. In contrast, under the same conditions but with 
$\mathcal {A}$ increased to 1 % as seen in figure 8(b), the turbulence puff now stops moving and stabilizes in about half a flow-through time, resulting in a significantly earlier reattachment length of 
$L_R/h \approx 67$. This behaviour with a steady state achieved after approximately one flow-through time was observed for all turbulent and laminar simulations presented in this paper.
Evolution of the centreline instantaneous streamwise velocity 
$u_x$ with inlet perturbation amplitudes of (a) 0.01 % and (b) 1 %, maximum 
${\rm CFL} = 0.5$ and grid G4 at 
$Re = 1500$. Flow from left to right (in regards to 
$x$) and from bottom to top (in regards to 
$t$).
3.2. Grid resolution
Figure 9 shows time-averaged streamwise velocity 
$u_x$ results for different grids along the centreline downstream of the expansion, for simulations without inlet perturbations, at 
$Re = 2500$. Similar to figure 8, it is possible to identify where the recirculation region ends and where the turbulent puff is positioned through the sharp drop in the centreline velocity as the flow moves downstream. The finer grids G4 to G1 all resulted in the same fully laminar state in these simulations. When testing coarse grids G5 and G6, transition to turbulence was observed, with the coarser grid G6 resulting in earlier transition to turbulence than G5. The main source of the perturbation that triggers and sustains turbulence here appears to be the numerical error from using coarser grids with larger spatial discretization errors in the solution of the equations, resulting in simulations with clear grid dependence when comparing G6 and G5 to the others.
Grid test centreline time-averaged streamwise velocity 
$\langle u_x \rangle$ results without inlet perturbations, with maximum 
${\rm CFL} = 0.5$ at 
$Re= 2500$. Flow from left to right.
On the other hand, when adding significant inlet perturbations, as seen in figure 10, all grids result in the formation of a turbulence puff at a similar position downstream of the expansion, and the grid dependence of the results is greatly reduced compared to without perturbations or with a fully laminar inlet as in figure 9. Also note that between figures 9 and 10, the results of G6 without or with inlet perturbations are effectively the same, indicating that a threshold of perturbation intensity was reached beyond which additional perturbations do not make a significant difference in the reattachment length results.
Grid test centreline time-averaged streamwise velocity 
$\langle u_x \rangle$ results with inlet perturbations at 
$\mathcal {A}=0.7\,\%$, with maximum 
${\rm CFL} = 0.5$ at 
$Re=2500$. The laminar curve indicates the laminar solution obtained when the simulation is unperturbed.
Further discussion of the grid resolution and a comparison to the Kolmogorov length scale is available in § 3.4.
3.3. Time-stepping resolution
A similar behaviour to the grid resolution is observed when testing different time steps. Overall, figures 11 and 12 show that the coarser (larger) the time step, the closer the location of the turbulence puff to the expansion or the earlier the transition from the laminar flow to turbulence. Under these conditions, for time steps with a maximum CFL of 0.5 and higher, a threshold seems to be reached, as the results with a maximum CFL of 0.8 and 1.0 are very similar to those with a maximum CFL of 0.5. On the other hand, these results also indicate that the time step may need to be smaller than CFL 0.1 in order to obtain DNS that are independent of time-step sizes, as the results are still changing around CFL 0.1.
Time-stepping results without inlet perturbations for centreline (a) time-averaged streamwise velocity and (b) Reynolds streamwise normal stress, with grid G6 at 
$Re= 2500$.
Time-stepping centreline time-averaged streamwise velocity results with inlet perturbations of 
$\mathcal {A}=0.7\,\%$, and grid G4 at 
$Re= 2500$.
Additionally, figure 11(b) shows the Reynolds streamwise normal stress results for these under-resolved simulations, which is an indication of the intensity of the streamwise velocity fluctuations along the centreline of the pipe. The peak value of the fluctuations progressively decreases as the time-step size is decreased. This indicates reduced turbulence intensity or a less disturbed simulation, and agrees with the notion that numerical discretization errors are perturbing the simulation. This reduction in velocity fluctuation peaks could also be explained by considering that smaller time steps will more accurately capture the behaviour of the smaller flow scales, which are responsible for turbulence dissipation. In particular, the results for CFL 0.05 depict a less disturbed simulation around 
$x/h=50$, but the growth of the fluctuations is slightly faster than for CFL 0.1 as seen by the inclination of the curve around 
$x/h=100$, though the peak value of fluctuations for CFL 0.05 is still smaller than for CFL 0.1. The change from the simulation with CFL 0.1 to the one with CFL 0.05 is much smaller than that from CFL 0.2 to CFL 0.1, even though the time-step size was halved in both cases, indicating that this may be a case of oscillatory convergence. In other words, at this point the solution seems to be converging with some oscillations to a time-step-independent solution.
The time-step sizes used here were significantly smaller than the Kolmogorov time scales 
$(\nu / \varepsilon )^{1/2}$ in these simulations. For example, with grid G6 and CFL 0.8 at 
$Re = 2500$ (figure 11), the simulation time-step size was about five times smaller than the smallest time scales calculated for each control volume in the grid. In comparison, the study of Choi & Moin (Reference Choi and Moin1994) investigated the influence of the time-step size on plane channel flow at a higher 
$Re$ (based on the channel half-width) of 4200, with time steps corresponding to CFL numbers up to 5 and closer to the size of the Kolmogorov time scale. The authors observed that at such large time steps, the simulation would result in a laminar solution. On the other hand, when using time steps of the order of CFL 0.5 to 3, a larger time step would result in larger streamwise velocity and normal vorticity fluctuations, qualitatively agreeing with the behaviour observed in the present study.
Lastly, as figure 12 shows, when adding inlet perturbations the time-stepping dependence of the results is greatly reduced, which is very similar to what was seen for grid dependence in figure 10. In general, there is surprisingly little difference between the streamwise velocity time-averaged results of these simulations with perturbations and different grids or time steps (figures 10 and 12, respectively), considering that the spatial resolution changed by more than 20 times from G1 to G6 and the temporal resolution changed by around 10 times from maximum CFL 0.1 to 1.0.
3.4. Turbulence kinetic energy budget
As mentioned previously, numerical discretization errors were observed to influence the simulation in a similar manner to random perturbations. Additional understanding of the behaviour observed in figures 9 to 12, with earlier transition to turbulence obtained when using a coarser discretization in time and space, comes from consideration of the role of the smaller-scale motions or eddies in the turbulence energy cascade. By using coarser grids and larger time steps, the smaller-scale motions are not adequately captured. Knowing that the smaller scales are associated with turbulence dissipation, it follows that, without the smaller scales to counteract the growth of turbulence, the simulation more readily grows into turbulence and maintains the turbulent state.
In order to substantiate this discussion, the cross-section-averaged turbulence kinetic energy budget was evaluated along the length of the pipe downstream of the expansion. Considering the Reynolds decomposition as 
$\boldsymbol {u} = \langle \boldsymbol {u} \rangle + \boldsymbol {u}'$, with the time-averaged velocity 
$\langle \boldsymbol {u} \rangle$ and the velocity fluctuation 
$\boldsymbol {u}'$, the equation of the mean turbulence kinetic energy budget 
$k = \frac {1}{2} \langle \boldsymbol {u}' \boldsymbol {\cdot } \boldsymbol {u}' \rangle$ can be written as (Pope Reference Pope2000, p. 125)
\begin{equation} \frac{\partial k}{\partial t} = \mathcal{P} - \langle \boldsymbol{u} \rangle \boldsymbol{\cdot} \boldsymbol{\nabla} k - \varepsilon- \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{T} , \end{equation}
where 
$\mathcal {P}$ is the production of turbulence kinetic energy calculated as
\begin{equation} \mathcal{P} ={-} \langle \boldsymbol{u}' \boldsymbol{u}' \rangle \boldsymbol{:} \boldsymbol{\nabla}\langle \boldsymbol{u} \rangle , \end{equation}
where 
$\boldsymbol {:}$ indicates the scalar product (or double dot product) operation.
The dissipation of turbulence kinetic energy 
$\varepsilon$ is calculated as
\begin{equation} \varepsilon = 2 \nu \langle \boldsymbol{{s}} \boldsymbol{:} \boldsymbol{{s}} \rangle , \end{equation}
where 
$\boldsymbol {{s}}$ is the fluctuating rate-of-strain tensor, given by
\begin{equation} \boldsymbol{{s}} = \tfrac{1}{2} ( \boldsymbol{\nabla} \boldsymbol{u}' + (\boldsymbol{\nabla}\boldsymbol{u}')^\mathrm{T} ) . \end{equation} The term 
$\langle \boldsymbol {u} \rangle \boldsymbol {\cdot } \boldsymbol {\nabla } k$ represents the mean-flow convective transport of 
$k$. Finally, the vector 
$\boldsymbol {T}$ in the diffusion or turbulence transport term is calculated as
\begin{equation} \boldsymbol{T} = \tfrac{1}{2} \langle ( \boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{u}' ) \boldsymbol{u}' \rangle + \langle p' \boldsymbol{u}' \rangle - 2 \nu \langle \boldsymbol{{s}} \boldsymbol{\cdot} \boldsymbol{u}' \rangle , \end{equation}
where 
$p'$ is the fluctuating kinematic pressure (i.e. 
$p' = p - \langle p \rangle$). To represent the imbalance or residual of the turbulence kinetic energy budget, a term represented by 
$E$ is subtracted from the right-hand side of (3.1) to ensure that 
$\partial k / \partial t$ is equal to zero.
A ratio between the control volume edge length 
$\varDelta$ and the Kolmogorov length scale 
$\eta$ was calculated for each control volume or cell in the grid, and will be represented as 
$\varDelta / \eta$. Each control volume's edge length is calculated by approximating the volume as a cube, and the Kolmogorov length scale for each volume is calculated as
\begin{equation} \eta = \left(\frac{\nu^3}{\varepsilon} \right)^{{1}/{4}}. \end{equation} Figure 13 shows the cross-section-averaged results for 
$Re = 800$ of the turbulence kinetic energy budget terms as a function of the location downstream of the expansion, with results non-dimensionalized and scaled by multiplying them by 
$10^{3} d/U^3$. In figure 13, laminar flow is coming from the left, and, after the turbulence puff, laminar flow leaves to the right. At the upstream front, coming from laminar flow, there is first a sharp increase in turbulence production. Similar to the turbulence slugs in a straight pipe evaluated by Wygnanski & Champagne (Reference Wygnanski and Champagne1973) and Song et al. (Reference Song, Barkley, Hof and Avila2017), there is a delay in the dissipation peak compared to the production peak. According to Wygnanski & Champagne (Reference Wygnanski and Champagne1973), this behaviour is explained by considering that only the smallest eddies contribute to the dissipation and therefore a finite time is required for the cascading process to properly take effect.
Turbulence kinetic energy budget results, with grid G4 and maximum 
${\rm CFL} = 0.5$ at 
$Re= 800$ with 
$\mathcal {A}=20\,\%$. Flow from left to right.
After the peak of production, as the flow approaches the downstream front of the puff, it moves back towards a laminar state with a decrease in both the turbulence production and the dissipation. In contrast, in the results of Song et al. (Reference Song, Barkley, Hof and Avila2017) for a slug at 
$Re = 5000$ (noting that the slug expands and moves with the flow whereas the puff here is fixed in size and position), there are secondary production and dissipation peaks near the downstream front of the slug. These downstream peaks are not seen at 
$Re = 2000$ in their slug results (referencing this behaviour as a weak front) nor in any of the results obtained in the present study for the puff in a pipe with an expansion for 
$Re \leqslant 2500$.
Still looking at figure 13, note that the term 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {T}$ has a very small contribution compared to the other terms, which is similar to what was observed for the turbulence slugs by Song et al. (Reference Song, Barkley, Hof and Avila2017). In other words, turbulence diffusion is minimal in these cases compared to its convective transport. The convective transport term 
$-\langle \boldsymbol {u} \rangle \boldsymbol {\cdot } \boldsymbol {\nabla } k$ also shows that convection moves turbulence away from the upstream front towards the downstream front of the turbulence puff.
Figure 14 presents an example of the turbulence kinetic energy budget calculated for a simulation with significantly reduced discretization resolution compared to figure 13, as the Reynolds number is increased from 800 to 2500 simultaneously with a grid significantly coarser. Note that the coarse simulation for figure 14 did not require the addition of inlet perturbations in order to produce turbulence, while the simulation in figure 13 did. Two features are immediately evident in the coarse simulation compared to the previous one: (i) the budget imbalance represented by 
$E$ is significantly increased, and (ii) the calculated resolved dissipation 
$\varepsilon$ is reduced compared to the production 
$\mathcal {P}$.
Coarse simulation results for the turbulence kinetic energy budget, with grid G6 and maximum 
${\rm CFL} = 0.8$ at 
$Re= 2500$ without inlet perturbations.
The reasons for the aforementioned behaviour are related: (i) the increased budget imbalance 
$E$ is due to the underestimated dissipation 
$\varepsilon$, and (ii) the dissipation is underestimated due to the coarse discretization with large time steps and large grid cells that cannot resolve the smaller scales, which are the ones responsible for turbulence dissipation. In this case, the numerical dissipation from the simulation algorithms seems to compensate for the failure to resolve the smallest scales, and contributes a sink equivalent to the budget imbalance 
$E$. In other words, 
$E$ comes from numerical or artificial dissipation from the solver itself. If numerical dissipation did not compensate for the lack of resolved dissipation, the simulation would probably diverge.
As a way to quantify the previous observations, the grid resolution was evaluated through the ratio between each control volume's edge length and the Kolmogorov length scale, 
$\varDelta / \eta$, see (3.6). It is generally recommended that the grid spacing be of the same order of magnitude as the Kolmogorov length scale for DNS, with Yeung & Pope (Reference Yeung and Pope1989) recommending 
$\varDelta / \eta$ up to 2.0 and Shin, Sandberg & Richardson (Reference Shin, Sandberg and Richardson2017) up to around 2.5. Here, the simulation for figure 13 resulted in 
$0.06 \leqslant \varDelta / \eta \leqslant 2.2$. In this simulation, an estimated 92 % of the dissipation was resolved (i.e. calculated through 
$\varepsilon$) and 8 % of the dissipation was associated with numerical dissipation (assuming and approximating that it is equal to 
$E$). In contrast, the coarse simulation in figure 14 resulted in 
$0.08 \leqslant \varDelta / \eta \leqslant 9.5$, with around 40 % of the dissipation coming from 
$\varepsilon$, and 60 % from 
$E$. The simulation for figure 13 was clearly of much higher resolution than the one for figure 14.
The approach of performing DNS with less than ideal resolution is not new. Komen et al. (Reference Komen, Camilo, Shams, Geurts and Koren2017), for example, using effectively the same OpenFOAM second-order-accurate numerical set-up that was used here and described previously, investigated fully developed turbulent channel flow. Of note, in a highly resolved simulation that the authors referenced as ‘quasi-DNS’, they also observed numerical dissipation responsible for 8 % of the total dissipation, which is the same value as reported here for the simulation shown in figure 13. For this simulation, Komen et al. (Reference Komen, Camilo, Shams, Geurts and Koren2017) report DNS-quality results in comparison to reference DNS databases. The authors also compared the usage of even coarser ‘under-resolved’ DNS with numerical dissipation estimated to account for 30 % to 60 % of the total dissipation, as an alternative to simulations with subgrid-scale modelling in the large-eddy simulation approach, and found the differences to be within a few percentage points. For related discussion on numerical dissipation and its effects on the simulation, interested readers are referred to Zhou et al. (Reference Zhou, Grinstein, Wachtor and Haines2014) and Cadieux, Sun & Domaradzki (Reference Cadieux, Sun and Domaradzki2017).
Figure 15 shows a simulation with intermediate resolution between those of figures 13 and 14 with 
$0.01 \leqslant \varDelta / \eta \leqslant 4.6$. This resulted in an estimated 70 % of the dissipation coming from 
$\varepsilon$, and 30 % from 
$E$, which also fits between the results of the other simulations. Comparing figures 14 and 15, which are at the same Reynolds numbers, more closely, the non-dimensionalized peak value of the turbulence production in figure 15 was half that of figure 14. This is despite the fact that figure 15 had significant perturbations added to the inlet, while figure 14 did not have any additional perturbations. This behaviour shows how much of a disturbance the numerical errors from the discretization can cause in the simulations, even triggering and sustaining turbulence by themselves, without the additional perturbations being required for turbulent flow.
Turbulence kinetic energy budget results, with grid G1 and maximum 
${\rm CFL} = 1.0$ at 
$Re= 2500$ with 
$\mathcal {A}=0.7\,\%$.
In summary, the main observations from this subsection are as follows: (i) Not resolving the smaller flow scales by using coarser discretization procedures results in reduced resolved dissipation 
$\varepsilon$ compared to the production 
$\mathcal {P}$. (ii) The simulation seemingly compensates for the reduced resolved dissipation with increased numerical dissipation. (iii) Finally, coarser discretization procedures also result in numerical errors that are a significant source of disturbance to the simulations, sufficient to trigger and sustain turbulence by themselves.
3.5. Inlet perturbations
When testing increased amplitudes of the arbitrary perturbations added to the inlet boundary of the simulation, the behaviour was similar to using coarser grids or time steps: larger inlet perturbation amplitudes result in significantly earlier transition to turbulence, as shown for example in figure 16. On the other hand, the influence of the amplitude of the inlet perturbations for a lower 
$Re$ of 800 is shown in figure 17, where a perturbation amplitude of 5 % was not large enough to trigger turbulence. The velocity fluctuations 
$\langle u'_x u'_x \rangle$ in figure 17(b) for 
$\mathcal {A}=5\,\%$ are not amplified downstream of the expansion and the flow remains fully laminar, in contrast with the results for larger perturbation amplitudes that were amplified and resulted in velocity fluctuation peaks and turbulence. For the simulation with 
$\mathcal {A}=10\,\%$, although it did exhibit a significant amplification of its velocity fluctuations, the streamwise velocity results in figure 17(a) are still closer to the fully laminar results obtained with 
$\mathcal {A}=5\,\%$. In contrast, the streamwise velocity fluctuations of the simulation with 
$\mathcal {A}=20\,\%$ are amplified by almost four orders of magnitude from 
$x/h \approx 25$ to 
$x/h \approx 70$. Lastly, throughout these simulations with different perturbation amplitudes 
$5\,\% \leqslant \mathcal {A} \leqslant 20\,\%$ at 
$Re = 800$, the time-averaged reattachment length did not change significantly from the laminar state around 
$L_R /h \approx 71.5\pm 3$ (see also figure 7).
Centreline time-averaged streamwise velocity profiles showing the influence of the inlet perturbation amplitude with grid G4 and maximum 
${\rm CFL} = 0.5$ at 
$Re= 2500$.
Centreline profiles for (a) time-averaged streamwise velocity and (b) Reynolds streamwise normal stress, with different inlet perturbation amplitudes, grid G4 and maximum 
${\rm CFL} = 0.5$ at 
$Re= 800$.
Figure 18 compares simulation results for different Reynolds numbers and perturbation amplitudes to results from the literature. The vertical axis shows the reattachment length, which represents the length of the recirculation region that is formed downstream of the expansion (see figure 1), obtained for different inlet Reynolds numbers shown on the horizontal axis. The experimental studies of Latornell & Pollard (Reference Latornell and Pollard1986) and Pak et al. (Reference Pak, Cho and Choi1990) used a sudden expansion ratio of 
$D/d =2$, while Back & Roschke (Reference Back and Roschke1972) used 
$D/d = 8/3$. In the numerical study of Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018) and in the present work, a ratio of 
$D/d=2$ was used. Pak et al. (Reference Pak, Cho and Choi1990) also performed experiments using 
$D/d = 8/3$, and the non-dimensional reattachment length results were effectively the same as in their 
$D/d =2$ experiments. The point where transition from laminar to turbulent flow starts for each study is the Reynolds number (flow rate) where the reattachment length peaks and begins to decrease. The steady laminar results for all studies, prior to the development of turbulence, closely agree with each other. The differences in the turbulent results can be explained by considering that the authors had different background disturbances or perturbation levels in their experiments (apart from other differences in the experimental apparatus), therefore exhibiting differences in the transition to turbulence. The numerical simulations of Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018) seem to be less subject to disturbances having a later transition to turbulence (higher 
$Re$ and higher reattachment lengths) compared to the other studies. Note that the simulation results of the present study using perturbations with 
$\mathcal {A} = 5\,\%$ closely agreed with the turbulent experimental results of Latornell & Pollard (Reference Latornell and Pollard1986). The concave curve shape formed by the turbulent results of the present study with constant 
$\mathcal {A}$ also agrees with the other studies.
Reattachment length (
$L_R$) results for different Reynolds numbers (
$Re$) and perturbation amplitudes (
$\mathcal {A}$), compared to literature results.
A more general interpretation and extrapolation of the results shown in figure 18 is schematized in figure 19. The shaded region indicates the multitude of possible turbulent flow states that can be produced depending on the perturbations affecting the flow. When sufficiently small or no perturbations are influencing the flow, it will result in a reattachment length given by the laminar flow line. It is possible that there is a lower boundary of Reynolds number below which flow is always laminar regardless of perturbation amplitude. Note that, for this 
$L_R$ versus 
$Re$ graph, it is also possible that the turbulent solutions region collapses towards the laminar result, with the turbulent simulations (provided with sufficiently small perturbations) resulting in essentially the same reattachment length as the laminar solution. This was seen, for example, in the behaviour of the simulations at 
$Re = 800$ (see the discussion for figures 7 and 17). As the aforementioned boundary, where turbulent states do not exist or result in the same reattachment length as the laminar ones, is likely to depend on the nature of the perturbation being applied, and the form and existence of such a boundary is not yet clear, it was not included in figure 19.
Suggested reattachment length chart with possible laminar and turbulent states.
Figure 20 presents the simulation results observed to be close to the critical point of transition from laminar to turbulent flow, when evaluating the effect of the inlet perturbation amplitude. A power-series fit of the results closer to the fully laminar limit or threshold, above which increasing the perturbation amplitude was observed to disrupt the laminar flow, was found to be approximately proportional to 
$Re^{-9.6}$. In contrast, Nguyen et al. (Reference Nguyen, Shadloo, Hadjadj, Lebon and Peixinho2019) tested by numerical simulations the amplitude of asymmetric vortex perturbations at the inlet boundary and observed a scaling with 
$Re^{-3}$ for 
$1000\leqslant Re\leqslant 2000$; Lebon et al. (Reference Lebon, Nguyen, Peixinho, Shadloo and Hadjadj2018a) tested by physical experiments a crosswise or transverse jet injection at the wall before the expansion for 
$400\leqslant Re\leqslant 1000$, and observed a scaling with 
$Re^{-2.3}$; Sanmiguel-Rojas & Mullin (Reference Sanmiguel-Rojas and Mullin2012) numerically tested ‘tilt perturbations’ at the inlet, consisting of small perturbations added only to the transverse component of the velocity vector with the other components unchanged, and observed a scaling with 
$Re^{-0.006}$ for 
$1300\leqslant Re\leqslant 2600$. Clearly, the 
$Re$ scaling of the critical perturbation amplitude seems to be significantly influenced by the characteristics of the perturbation being evaluated, as well as the range of 
$Re$ tested.
Behaviour of simulation results for different 
$Re$ and perturbation amplitudes close to the critical point of transition: ‘laminar’ indicates fully laminar flow with the absence of a turbulence puff; ‘turbulent’ indicates the existence of a stable turbulence puff; and ‘unstable’ indicates simulations that had neither a stable turbulence puff nor fully laminar flow.
A turbulent solution state was not found at 
$Re = 500$ when testing perturbation amplitudes up to 50 %, a result that qualitatively agrees with that of Lebon et al. (Reference Lebon, Nguyen, Peixinho, Shadloo and Hadjadj2018a), who estimated a lower bound of 
$Re$ around 400 below which only laminar flow would be observed. Peixinho & Mullin (Reference Peixinho and Mullin2007) also found a lower 
$Re$ bound for turbulent pipe flow without an expansion. It is possible that these lower bounds observed can result from (i) not testing perturbations with amplitudes high enough, or (ii) limitations due to the nature of the different types of perturbations tested. For reference, the power-law fit shown in figure 20 indicates that perturbation amplitudes 
$\mathcal {A} > 100\,\%$ would be required in order to observe turbulence for 
$Re \le 500$.
3.6. Inlet disturbance behaviour
The previously discussed figures 11(b) and 17(b) showed the behaviour, downstream from the expansion, of the disturbances through the centreline time-averaged velocity fluctuations for simulations without and with added inlet perturbations, respectively. In order to get a better understanding of how the disturbances behave before the expansion, i.e. in the smaller-diameter inlet pipe, a different set of figures is presented.
Figure 21 shows the results of a simulation that did not require added inlet perturbations in order to produce the turbulence puff. The flow disturbances, being evaluated as velocity fluctuations, more specifically through the magnitude of the time-averaged Reynolds stress tensor, are negligible at the inlet (
$x/h=-10$), of the order of 
$10^{-14}$, as seen in figure 21. This is expected, as this simulation did not have any added inlet perturbations, with a laminar profile specified at the inlet. The only disturbance at the inlet is numerical noise, likely to be due to truncation and rounding errors. However, as the flow moves through the smaller-diameter inlet pipe, from 
$x/h=-10$ to 
$-5$, there is a significant growth of these disturbances, reaching values around 
$10^{-4}$ at 
$x/h=-5$. Afterwards, there is no immediate and significant change in the magnitude of the velocity fluctuations, from upstream of the expansion, 
$x/h=-5$, to downstream, 
$x/h=5$. However, as the flow moves downstream further along the expanded region of the pipe, with the recirculation region along the walls and an adverse pressure gradient, these disturbances eventually grow to a point where turbulent flow is observed, with values of 
$\lvert \langle \boldsymbol {u}'\boldsymbol {u}' \rangle \rvert$ up to 
$10^{-1}$ at the reattachment length of 
$L_R = x/h = 93.2$.
Time-averaged inlet pipe velocity fluctuations for 
$Re = 2500$ and grid G5, without added inlet perturbations.
On the other hand, looking at a simulation that did require perturbations in order to produce turbulence in figure 22, a completely different behaviour is observed. Initially, at the inlet (
$x/h=-10$), the added inlet perturbations produce significant velocity fluctuations with values around 
$10^{-2}$. However, as the flow moves through the inlet pipe to 
$x/h=-10$, these decay to less than 
$10^{-4}$. It is interesting to note that, although there is no significant asymmetry in the fluctuations due to the perturbations at the inlet (
$x/h=-10$), there is a notable amount of asymmetry in the profiles closer to the expansion, as seen in 
$x/h=-5$. This behaviour was not observed in the results in figure 21. Then, similar to the previous result, as the flow moves downstream from the expansion, the disturbances grow, and in the region where the flow reattaches, velocity fluctuations reach values around 
$10^{-1}$, significantly surpassing the magnitude of the fluctuations from the perturbations at the inlet.
Time-averaged inlet pipe velocity fluctuations for 
$Re = 1500$, grid G4 and 
$\mathcal {A}=5\,\%$.
Figure 23 compares (a) a simulation with perturbations of 
$\mathcal {A}=5\,\%$, that were not enough to observe turbulent flow, to (b) a simulation with perturbations of 
$\mathcal {A}=20\,\%$, where the perturbations were enough to trigger turbulence downstream from the expansion. In this figure, it is immediately noticeable that reducing the inlet perturbation amplitude from (b) 
$\mathcal {A}=20\,\%$ to (a) 
$\mathcal {A}=5\,\%$ causes all velocity fluctuations from 
$x/h=-10$ to 
$x/h=5$ to be reduced by at least one order of magnitude. On the other hand, the shape of the curves from 
$x/h=-10$ to 
$x/h=5$ is maintained when the perturbations are reduced from 
$\mathcal {A}=20\,\%$ to 
$\mathcal {A}=5\,\%$, with the same asymmetries. The inlet perturbations of 
$\mathcal {A}=5\,\%$ resulted in velocity fluctuations close to 
$10^{-2}$ at the inlet, but only laminar flow was observed with fluctuations of less than 
$10^{-4}$ at 
$x/h=71.5$, which is close to where the flow reattaches. On the other hand, with inlet perturbations of 
$\mathcal {A}=20\,\%$ (figure 23b), fluctuation magnitudes closer to 
$10^{-1}$ are observed at the inlet, and downstream of the expansion this resulted in the formation of the turbulence puff with a similar magnitude of velocity fluctuations to the inlet, above 
$10^{-2}$.
Time-averaged inlet pipe velocity fluctuations for 
$Re = 800$, grid G4, and (a) 
$\mathcal {A}=5\,\%$ and (b) 
$\mathcal {A}=20\,\%$.
Lastly, looking at figures 21–23 as a whole, the simulations that produced turbulence (all except for figure 23a) had similar non-dimensional velocity fluctuation magnitudes close to 
$10^{-1}$ at the reattachment region (which is also where the turbulence puff is located). At these values, visible turbulence and chaotic flow are observed, whereas at fluctuation magnitudes below around 
$10^{-3}$ the flow was effectively still laminar, steady and organized (also see figures 7 and 17, which are under the same conditions as figure 23).
4. Limitations and suggestions for future work
In the present work, as mentioned previously, the length of the pipe after the expansion was limited to 
$L_o = 200 h$. In comparison, Moallemi & Brinkerhoff (Reference Moallemi and Brinkerhoff2018) used a geometry with 
$L_o = 100 h$ and Vittal Shenoy et al. (Reference Vittal Shenoy, Safdari Shadloo, Peixinho and Hadjadj2019) one with 
$L_o = 300 h$. Overall, in order to reduce the required computational resources, the fluid domain sizes were limited in both a spatial as well as a temporal sense. We believe these could be important for the simulations closer to the point of global transition to turbulence (e.g. unstable simulations in figure 20), where the behaviour is more strongly intermittent, though these were not the focus here. At least for the simulations that produce a stable turbulence puff or that result in a fully laminar profile, the geometry and time simulated were sufficient to observe fully developed steady-state behaviour (see for example figure 8 and its discussion).
The length of the smaller-diameter inlet pipe can play an important role in the behaviour of this system, but its influence was not evaluated here. A longer inlet pipe would cause perturbations coming from the inlet boundary to decay further before they are amplified by the sudden expansion (provided they are not enough to trigger turbulence in the inlet pipe itself). In other words, a longer inlet pipe would further diminish the influence of perturbations at the inlet. It follows that, if the length of the inlet pipe was longer, higher inlet perturbations would be necessary in order to compensate for it and achieve the same effect, and vice versa if it was shorter.
The simulations with different Reynolds numbers were performed using fixed grids. Effectively, this is expected to result in higher spatial discretization errors for simulations with higher Reynolds numbers when compared to simulations with the same grid but lower Reynolds numbers. Here, the grid convergence was evaluated mainly for the highest Reynolds number used of 2500 (see figures 9 and 10). Ideally, a separate grid convergence study would be performed for each different Reynolds number being simulated, in order to better evaluate the influence of grid refinement and numerical perturbations coming from the spatial discretization errors. Additionally, in terms of grid dependence, only the global resolution or refinement of the grid was evaluated here (i.e. size of the grid cells). However, changes in the construction, such as in different cell polyhedron types, cell orientation, aspect ratios and/or refinement regions, would probably also be relevant. Lastly, simulation resolution was also evaluated using the calculated Kolmogorov length scales, as discussed in § 3.4.
The sources of perturbation used here to trigger turbulence were not expected to represent realistic behaviour, particularly considering that numerical error itself was the only source of perturbation in some of the turbulent simulations. Even though the results of the simulations overall agree with experimental and numerical results from the literature (see for example figure 18), additional simulations using more realistic inlet perturbations, extracted from physical measurements or from precursor DNS, could complement these results. In fact, the review of Mullin (Reference Mullin2011) noted that drawing a connection between the disturbances in experimental and numerical investigations was an open problem in the literature, and to the best of our knowledge this is still an open problem today.
The mechanisms behind how numerical errors can grow as disturbances and result in perturbations to the flow are unclear. It is possible that different perturbations lead to turbulence through different mechanisms. In addition, a technique to quantify numerical errors in terms of the equivalent perturbation amplitude that they cause could be useful in this discussion.
5. Conclusions
The most relevant differences in the turbulent behaviour of simulations and real-world physical experiments of pipes with a sudden expansion are associated with the fact that the simulations typically do not take into account sufficient and equivalent sources of perturbation in comparison to physical experiments. Physical experiments are susceptible to various imperfections such as wall scratches or defects that can cause perturbations and trigger transition to turbulence. In contrast, numerical simulations typically consider perfect or ideal conditions with an absence of imperfections, apart from intrinsic numerical errors. In the present paper, sources of numerical perturbation, including added inlet perturbations as well as numerical discretization errors, were evaluated.
Observations include that, at lower Reynolds numbers, such as 
$Re \leqslant 800$, turbulent states may not produce significantly different reattachment lengths compared to the laminar solution, even though the recirculation region changes significantly and becomes thicker before ending. At higher Reynolds numbers, such as 
$Re \geqslant 1000$, the recirculation region ends and the flow reattaches significantly earlier in the turbulent states. Additionally, the formation of a large helical vortex within the laminar recirculation region downstream of the expansion was observed in the turbulent solutions.
Four main mechanisms were identified to interact with each other and influence the flow in the simulation transitioning between laminar and turbulent states: (i) the smaller eddies resolved by the simulation; (ii) numerical dissipation from the simulation algorithms compensating for missing dissipation from unresolved smaller scales; (iii) spatial and temporal discretization numerical errors that disturb the simulations; and (iv) perturbations intentionally added to the simulations, with an example here of inlet boundary perturbations.
In summary, mechanisms (i) and (ii) seem to work together to provide a sink to the turbulence kinetic energy produced by the larger resolved scales, and mechanisms (iii) and (iv) work together to perturb the simulation and trigger the production of turbulence kinetic energy by the simulation. The influence of using coarser spatial and temporal discretization procedures, which cause larger numerical errors (iii), was found to be similar to the influence of adding or increasing the amplitude of inlet perturbations (iv). It is expected that mechanisms (i), (ii) and (iii) will always be present in any turbulent CFD simulation to some extent, regardless of how refined the resolution or methods of the simulation are.
Overall, despite the arbitrary inlet perturbations and the numerical errors not being expected to represent realistic sources of perturbation, the results are remarkably consistent between themselves and with results from the literature. They also reinforce the notion that the larger the Reynolds number, the more sensitive the simulation is to the amplification of perturbations and numerical errors that can result in turbulence. Perhaps counter-intuitively, this also implies that, in order to observe a laminar solution at higher Reynolds numbers, the simulation needs better numerics such as finer grids or smaller time steps, otherwise numerical errors by themselves can cause perturbations significant enough to be amplified and trigger turbulence in the flow, regardless of inlet or other sources of perturbations. On the other hand, if using fine grids and time steps to minimize numerical discretization errors, the addition of numerical perturbations of sufficient amplitude is necessary in order to trigger turbulence. Lastly, this demonstrates how reporting coarse simulations together with refined simulations can give insights into the simulation due to the differences in numerical disturbances and flow scales resolved.
Finally, turbulence that seems to arise naturally in numerical simulations has probably been triggered by the amplification of disturbances caused by numerical error. Perfectly unperturbed simulations without any numerical errors or disturbances would result in fully laminar solutions, regardless of Reynolds number. The main sources of numerical disturbances identified in the present study were from spatial and temporal discretization errors. Similar to how perfectly undisturbed physical experiments are not feasible, numerical errors in the simulations can be minimized, but it is unlikely that they can be completely eliminated. In other words, numerical errors are expected always to play a part in the transition to turbulence of numerical simulations.
Funding
The authors are grateful for financial support of the present study via scholarships and research funding, including the University of Saskatchewan Dean's Scholarship and the Devolved Scholarship from the Department of Mechanical Engineering, and research funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada. The authors also thank Research Computing at the University of Saskatchewan for the computational resources.
Declaration of interests
The authors report no conflict of interest.
Open access
