Features of Transition to Turbulence in Sudden Expansion Pipe Flows
Minh Quan Nguyen Benoit Lebon Mostafa Safdari Shadloo Abdellah Hadjadj Jorge Peixinho







UMR 6614-CORIA, Normandie University, CNRS-University and INSA of Rouen, 76000 Rouen, France


Laboratoire Ondes et Milieux Complexes, CNRS and Université Le Havre Normandie, 76600 Le Havre, France


UMR 6614-CORIA, Normandie University, CNRS-University and INSA of Rouen, 76000 Rouen, France


Laboratoire Ondes et Milieux Complexes, CNRS and Université Le Havre Normandie, 76600 Le Havre, France


UMR 6614-CORIA, Normandie University, CNRS-University and INSA of Rouen, 76000 Rouen, France


The complex flow features resulting from the laminar-turbulent transition (LTT) in a sudden expansion pipe flow, with expansion ratio of 1:2 subjected to an inlet vortex perturbation is investigated by means of direct numerical simulations (DNS). It is shown that the threshold for LTT  described by a power law scaling with -3 exponent that links the perturbation intensity to the subcritical transitional Reynolds number. Additionally, a new type of instability is found within a narrow range of flow parameters. This instability  originates from the region of intense shear rate which is a result of the flow symmetry breakdown. Unlike the fast transition, usually reported in the literature, the new instability emerges gradually from a laminar state and appears to be chaotic and strongly unsteady. Additionally, the simulations show a hysteresis mode transition due to the reestablishment of the recirculation zone in a certain range of Reynolds numbers. The latter depends on (i) the initial and final quasi-steady states, (ii) the observation time and (iii) the number of intermediate steps taken when increasing and decreasing the Reynolds number.








Features of Transition to Turbulence in Sudden Expansion Pipe Flows
     Abdellah Hadjadj
    December 30, 2023
===================================================================



§ INTRODUCTION


The flow through an axisymmetric sudden expansion in a circular pipe is a basic configuration, which occurs in many industrial applications, such as heat exchanger, mixing chamber, combustion chamber, etc.. This basic geometry is also used as a building block to model more complex flows such as those occurring in arterial stenoses  <cit.>, pistons <cit.>, and transportation pipes <cit.>, among others. In these applications, the capacity of predicting when the flow will become turbulent is crucial. In the literature, there are many efforts in theoretical analysis <cit.>, experimental explorations <cit.> and numerical simulations <cit.> focusing on this problem, or a similar geometries such as the planar expansion <cit.> and the abrupt contraction and expansion <cit.>. 
Sudden expansion with various expansion ratio in investigated by <cit.>. More recently, <cit.> found a new mechanism for periodic bursting of the recirculation region in the flow in a circular pipe with the expansion ratio of 1:2. Yet, a consensus about the sequence of event in the transition from laminar to turbulence seems to be relatively well-established, but the exact value of critical Reynolds number is still not firmly determined, and the different transition scenarios are not yet fully elucidated.

The flow is mainly controlled by the inlet Reynolds number, Re=Ud/ν, where U is the mean velocity at the inlet,  d is the inlet diameter and ν is the fluid kinematic viscosity. In the laminar state, the flow is axisymmetric. As Re increases, the flow starts to break its symmetrical properties but remains steady until Re = 1139± 10 as reported in the experiments by Mullin et al. <cit.>. Noted that the onset of symmetry breakdown is at much lower value for the case of channel sudden expansion, as reported to be Re≈40 by Fearn et al. <cit.> and Re=216 by Drikakis <cit.>. In the sudden circular pipe expansion flows, intermittent bursts were reported  <cit.> to appear at Re=635 and Re≈ 1453± 41. Then the flow enters to an oscillating state. For instances, this state found at Re≈ 1567 ± 16 <cit.>, 1500<Re<1700 <cit.>, and Re=750 <cit.> among others. Finally, the oscillating pattern breaks into a chaotic movement at even higher Reynolds numbers.

More recently, a global stability analysis was performed by Sanmiguel-Rojas et al. <cit.>. They showed that the flow can remain symmetry up to Re≈ 3273, which is much larger than the value found experimentally. Subsequent simulations by  Cliffe et al. <cit.> also indicated that the steady supercritical bifurcation point lies at even higher Reynolds numbers, i.e. Re≈ 5000. These discrepancies among the reported results suggest that there might be some missing underlying parameters. 

A detailed study of transient growth stability was performed by  Cantwell et al. <cit.>. They illustrated that the flow is linearly stable up to Re = 1400, however, even in the subcritical regime, the sudden expansion can amplify the energy of the perturbation up to 6 order of magnitude in the inlet and then decay. This difference may be explained by the fact that in experimental studies, as warned by <cit.>, the imperfection of apparatus has a strong impact on the measured critical Reynolds number. Therefore, the values of critical Reynolds number seem to be dependant both on the perturbation nature and its amplitude. In this case, a numerical simulation with a well-defined finite amplitude perturbation is required to better understand the underlying physics. 

The first direct numerical simulation (DNS) of finite amplitude perturbation in sudden expansion flow was recently performed by Sanmiguel-Rojas & Mullin <cit.> where many interesting results are reported. They showed that for a given Reynolds number there exists a region where the flow in laminar state can be forced to enter turbulence state with a finite perturbation. The minimal value of perturbation required to trigger the turbulence, scales accurately with Re^-0.006. This result has an important role in the passive flow control, however, it arises many questions. For instances, is this result universal? Will another kind of perturbation have the same behaviour? Is the scaling law still valid? The question is crucial especially by knowing that the previous perturbation scheme used by Sanmiguel-Rojas & Mullin <cit.> is simply a tilt disturbance in the velocity field. Therefore, it does not respect the no-slip boundary condition at the inlet wall. Another interesting result is: when they increased and then decreased the Reynolds number, a hysteresis behaviour was found for the Reynolds numbers ranging from 1450 to 1850. This result has strong implications on the nature of transition. This important point deserves more attention, since in the original work, the process of variation of Reynolds number as well as the physical time of the reported flow state are not specified. One may ask, will the hysteresis not occur if the Reynolds number vary in quasi-static manner? Or will the results change if the observation time is different? In the present study, the authors propose a similar study, but with vortex disturbance, to revisit the power law and the hysteresis behaviors. 



§ NUMERICAL SET-UP

The present work focuses on a pipe flow with a sudden expansion with the expansion ratio of 1:2. The fluid flow system is solved using Nek5000 <cit.>, a well-validated high-order spectral element code for transitional and fully turbulent flows <cit.>. The governing equations are mass and momentum conservations in an isothermal incompressible limit:


    ∇·𝐮=0



    ∂/∂ t𝐮+𝐮·∇𝐮=-∇ p+νΔ𝐮

Where 𝐮 is the velocity field, p is the pressure, and t is the physical time. The density of flow is set to unity for simplicity.

For the expansion pipe test case, the computational domain is axisymmetric (see figure <ref>). The region upstream of the expansion is called inlet and has the diameter of d and the length of 5d . The downstream is the region where we are focusing our analyse on. As illustrated in figure <ref>, this region has the diameter of D=2d  and the length of L=150d. The whole domain contains 62 300 spectral elements where each element consists of P^3 Gausse-Legendre-Lobatto (GLL) points, P being the polynomial order.  The element with P=5, same as <cit.>, with total number of 7.9 million calculating points is used in the whole domain for the reported results. However, to confirm our finding, additional simulations are performed with P=6 (13.5 million GLL points). A classical set up for the inlet velocity, located at z=-5d, is the Hagen-Poiseuille profile, which should satisfy the non slip condition at inlet wall: 𝐮_inlet^o(r=d/2,z=-5d)=0 when r=√(x^2+y^2)=d/2 : 







    𝐮_inlet^o(x,y,z=-5d)=2[1-4(x^2+y^2)/d^2]  U𝐞_𝐳,

where (𝐞_𝐱,𝐞_𝐲,𝐞_𝐳) are the set of 3 unit vectors of Cartesian base, respectively. The length unit  d and the time unit s are chosen such that the mean inlet velocity U reproduces unity. 

In order to trigger turbulence in subcritical Reynolds numbers, a vortex perturbation is added to the inlet parabolic profile. This modifies the expression of velocity inlet as:

    𝐮_inlet(x,y,z=-5d)=𝐮_inlet^o+𝒜Ω( [ - y;   x;   0 ])=( [               -𝒜Ω y;                𝒜Ω x; 2(1-4(x^2+y^2_r^2)) ]),







where 𝒜 and Ω are the amplitude and the intensity of the vortex perturbation, respectively. The value of Ω is controlled by its radius, R_Ω=d/4, and its center (x_Ω,y_Ω). Let r_Ω be the distance from a point (x,y) to the center of the vortex: r_Ω=√((x-x_Ω)^2+(y-y_Ω)^2), then the expression of Ω(x,y) is given by:

    Ω(x,y)= {[               1      r_Ω≤ 0.125;       2(1-4r_Ω) 0.25≥ r_Ω>0.125;               0        r_Ω>0.25 ]. .

The parameter of Ω(x,y)  is fixed with x_Ω=-0.25d and y_Ω=0, 
 such as the inlet profile respect the incompressible, non slip and non penetrate boundary condition. Indeed, the incompressibility is verified by: 

    ∇·u_inlet= {[       -∂(𝒜 y)/∂ x+∂(𝒜 x)/∂ y +  2∂ (1-4(x^2+y^2))/∂ z = 0                                                r_Ω≤ 0.125; 2∂(𝒜 yr_Ω)/∂ x-2∂(𝒜 xr_Ω)/∂ y +  2∂ (1-4(x^2+y^2))/∂ z =0                                           0.25≥ r_Ω>0.125;                                 2∂ (1-4(x^2+y^2))/∂ z = 0                                                  r_Ω>0.25 ]. .

The radius and the position of the vortex is chosen such that Ω vanishes at the boundary. For the closest point of the boundary to the center of the vortex: x=-0.5d and y=0, we have r_Ω=0.25d, thus Ω(x,y)=0. By injecting Ω=0 into equation <ref>, one can easily verify that 𝐮_inlet vanishes at the boundary as well. 





 Preliminary tests with the vortex disturbance on the pipe centreline indicate the required amplitude to initiate turbulence was very large <cit.>. It corresponds to the limit case of the rotating Hagen-Poiseuille flow discharging into a sudden expansion <cit.>. By positioning the vortex at (x_Ω,y_Ω)=(-0.25d,0), the vortex breaks the axial symmetry of the flow and naturally deflects the recirculation region. 
The distribution of Ω(x,y) can be visualised in figure <ref>(a). 







§ DISTINCTION BETWEEN DIFFERENT STATES OF INSTABILITIES


Laminar or turbulent nature of a flow can be clearly distinguished from the time evolution of the drag coefficient:
   
    C_z(t)=ν/4π d^2 L U^2∫_z=0^L∫_θ=0^2π[∂u_z/∂ r]_r=dr    dθ  dz,

where d and L are, respectively, the radius and the length of the downstream pipe. Here r=√(x^2+y^2) and θ=arctan(y/x) are the positions in the radial and azimuthal direction in the Cylindrical coordinate system, respectively.

Table <ref> presents a summary of the selected simulations along with their Reynolds numbers, perturbation intensities, initial conditions and observed states. While the laminar state is noted as LS, there exist 2 types of unsteady states with distinguished amplitude levels and positioning of turbulence puffs. Therefore, they will be called by a more general name, i.e. unsteady state (US). There are two fundamental differences between LS and US when observed directly from C_z(t). Firstly, LS is steady, therefore, the time variation of the skin-friction coefficient should be negligible, i.e. dC_z/dt≈ 0. This is not the case for US. Secondly, LS has a significantly smaller drag coefficient with regards to its US counterpart, i.e. C_z,LS<<C_z,US. 


To better observe the unsteady patterns, the velocity fluctuations are extracted in the pipe centerline by subtracting the instantaneous field u_z from a reference field u_zo in the same spatial location, i.e. 

    u'_z(0,0,z,t)=u_z(0,0,z,t)-u_zo(0,0,z),

and plotted in a space-time diagram. The reference value u_zo is the “closest" available laminar state for a given sets of parameters. In other words, u_zo is chosen as the initial condition or the final state of the simulation, if it starts or ends with a laminar state. However, if the simulation starts and ends with turbulence states, a laminar state from the closest Re and/or 𝒜 is chosen as the reference. 
The process of constructing such diagram is as the following: at a given instance, the value of interested quantities, here is u_z, of 1500 points in the centreline in Z direction is recorded at each 100 time-step. The final result is a 2D function of space and time with the resolutions of 0.1d and 0.1s), respectively. 
This diagram can provide complementary information to the plot of C_z, such as spatial distribution and local intensity of unsteady pattern. 

Figure <ref> presents the evolution of drag coefficient C_z as a function of time for four different perturbation amplitudes 𝒜 at Re=2000. The background colour in subfigures shows the space-time diagram. It is noted that the intensity of colour in the diagrams increases with the deceleration of the streamwise velocity. There are two main mechanisms that cause this deceleration. The first one is when the perturbation breaks the flow symmetry down. As a result, the centroid maximum value in the velocity profile is moved away from the centerline. This causes the streamwise velocity in the center of the pipe seemed to be decreased. If the flow is steady, but the symmetry is broken, a smooth region of blue colour will be observed. The second mechanism is when the flow becomes turbulent. As a consequence, the instantaneous streamwise velocity will fluctuate with strong amplitude. For the sake of simplicity, only deceleration is shown in the colormap. Additionally, the mean velocity profile is flattened when the flow becomes turbulent, which also results in a deceleration in the pipe centerline. Hence, the beginning of the turbulence has a clear signature in space-time diagram. This region is found with more intense blue colour, and a noisy interface in space-time diagram. Furthermore, it is obvious from this figure that the emergence of unsteady pattern is closely linked to the rise of drag coefficient C_z. Note that at the beginning of each simulation, if there is any change in the parameters comparing to the previous simulation, i.e. the simulation in which the last instance is used as initial condition for the new simulation, the flow might develop an instability. This instability gets carried downstream and decays. Since the cause of this instability is not physical, it will not be considered here. 




Several unsteady behaviours can be observed in the chosen range of parameters.
For all the cases, the flow exhibit a recirculation zone, which extent depend on Re. The recirculation zone is delimited by a reattachement point, the process of determination of this point is simple in laminar case: we find the first point near the wall where the velocity change the sign. For the turbulence cases, an average are over a brief periode of time is needed to determinate the point.  
For the case corresponding to Re=2000 and 𝒜=0.5 (see figure <ref>d), the unsteady behaviour starts almost immediately at the beginning of the simulation. The generation and sustention of the instability at the trailing edge seems to be linked with the breaking of the recirculation zone. The flow suddenly loses its smoothness and finds a new reattachment point, with the position right before the trailing edge. At 𝒜=0.5, this position is fixed over time. As the perturbation amplitude decreases, the flow has a tendency to restore the recirculation zone, and pushes the trailing edge further downstream. On the other hand, the unsteady pattern has tendency to profilate upstream and downstream, therefore try to push trailing edge upstream. The competition between these two mechanisms might lead to an unsteady position of the trailing edge such as that in the case of Re=2000 and 𝒜=0.2 (see figure <ref>(c)). The recirculation can be observed clearly with the contour plot of velocity in streamwise direction. 
As shown in figure <ref>(c), the unsteady pattern invades and pushes the reattachment point upstream compare to the laminar configuration (i.e. figure <ref>(a)). This type of unsteady behaviour will be called the primary unsteady-state pattern, hereafter labelled by US1. Since US1 breaks the recirculation region significantly and gains energy from this process, the transition from LS to US1 is irreversible and a hysteresis behaviour might be expected for this state.











For lower values of 𝒜, i.e. 𝒜≤ 0.09 (see figures <ref>(a) and <ref>(b)), the perturbation is weak and not able to spread out over the recirculation zone. Therefore, the unsteady pattern can not sustain or regenerate anymore; it gets carried downstream and the flow becomes globally laminar. However, due to the exerted perturbation at the inlet, a very small perturbation pattern with the bounded (i.e. small and quasi constant) amplitude still exists in the laminar state, which is called weak unsteady pattern (WUP). 
Depending on the values of 𝒜 and Re, WUP can stay bounded with a small amplitude (for instance, from t≈750s to t≈ 1000(s) in figure <ref>(b)), and appears globally as a LS. The flow could experience a transient growth and become a large scale unsteady state (for instance, from t≈1000 in figure <ref>(b)). The resulted state is clearly different from US1 described earlier and will be called as the second unsteady state pattern, hereafter denoted by US2. Comparing to US1, the unsteady pattern of US2 can evolved from an LS, the recirculation zone is mostly preserved (see figure <ref>(b)), and the overall amplitude of US2 is smaller than its counterpart in US1. An intermittency behaviour can be observed for the case of 𝒜≈ 0.0782 (see figure <ref>(a)), where the unsteady pattern of US2 emerged but decayed after. This pattern will be analysed in more details in the following sections. 




Other values of Reynolds numbers also show a similar behaviour. The observed states of the flow namely, LS, US1 and US2 are reported in Table <ref>. Note that while the state US1 can be firmly determined, the boundary between US2 and LS depend on how long we observe the system, since US2 can emerge from LS at very late in the simulation. Therefore, the state reported as LS could be US2 if the observation time is larger, and the boundary can not be firmly determined. Thus, the laminar state in Table <ref> are reported at least after t=1000(s). Larger observation time are used as Re approaching the critical value.






In the figure <ref>, the boundaries of LS and US2 are shown. The first observation is the border of transition follows a power law of 𝒜∼ Re^-3, which is much steeper than the value of -0.006 that is reported by Sanmiguel-Rojas & Mullin<cit.>. Meanwhile, the figure provided by the authors shows a completly different powerlaw, which estimated around -11.2. This important results yet confusing deserve more attendtion. Moreover, the perturbation method used by <cit.> is a tilt disturbance, i.e. shifting the inlet velocity profile upward (in X direction). This creates a discontinuity in the inlet and at the wall, which causes an intense shear in upstream flow. The effects of the shear generated could depend on the mesh resolution as well as on how the solver interpolate the discontinuity; but not only the physical parameters. A recent experimental study from <cit.> show a powerlaw of Re^-2.3, which is much closer to our results. 

The second observation is that the band of parameters that exhibits US2 is very narrow, and it situates in the small border between turbulence state US1 and the laminar state LS. Therefore, it plays an important role in understanding of the transition mechanism, especially when the Re and/or 𝒜 parameters evolves gradually. The state US2 is new compared to the state US1 which is well documented, so in the following section, this state will be analysed in more details. 














































































































































§ TRANSIENT GROWTH OF WEAK UNSTEADY PATTERN


In this section, the WUP state and the transition from LS to US2  will be further investigated. The smoothness of this transition make it particularly interesting to study, since it usually happens too fast to be captured.

Figure <ref>(a) presents the flow's space-time diagram for Re=1360 and 𝒜=0.2, corresponding to case 5a in Table <ref>.

This case is selected because, it contains unsteady regions and it is the closest case to an LS scenario. The flow stays for a long period in LS, and the transient growth arises around t≈ 2050.

The time evolution of the system can be divided into three phases. The initial phase when the flow adapts to the new parameters.  Before t<600, an un-physical instability develops and decays downstream. After the initial phase, the flow experiences a long lasting laminar state, for 600 <t<2000, with a bounded unsteady pattern WUP. Then, at t>2000 the flow  shows, for this particular case, an intermittent pattern between laminar state (LS) and turbulent state (US2). Note that For slightly higher Re, the turbulent state US2 is self sustained,  the intermittent pattern ceases.

 





To take a closer look at WUP, two time snapshots (figure <ref>b) of u'_z along the streamwise direction are taken at t=1000 (LS) and at t=2500(US2). The flow contains  WUP even in the laminar stage. Additionally, three distinct space regions, which two of them are identical in both LS and US2, can be seen.

The first zone is a steady region for 0<z<20d, where the flow is completely laminar and time independent. 
The second zone extent from the first zone to the reattachment point, approximately at 20<z<50d, where WUP show a clear wavy pattern in space with the wave length of λ_WUP≈ 2.34d. 



Figure <ref> presents two time signals of streamwise velocity fluctuations recorded respectively at z=40d and z=80d over the pipe centerline (i.e. x=y=0). The selected time signals in this zone (see figures <ref>(a) and <ref>(c)) have a similar noise's amplitude u'_WUP/U≈0.1 and frequency f_WUP≈0.2 Hz. This witnesses the similarities in the second zone for WUP in both LS and US2.

In the third and last zone, i.e. the region after the reattachment point, WUP transit into larger scale dynamics. In this region, the low frequency disturbances start to merge (see also figure <ref>(b) purple line), with larger wave length in the case of LS. If the parameters Re and δ are sufficiently high, the unsteady pattern in third zone can be amplified into a strongly time dependent and non linear pattern US2, see also figure <ref>(d). 

As mentioned earlier, the second and the third zones are situated before and after the reattachment point. This point, beside the mathematical definition, can also be visually determined from the contour plot of the streamwise velocity at level v_z=0 (see figure <ref>). The boundary between the second and the third zones is not sharp because the WUP shows organized structures which slowly decays as they crossed the border. This position indicates an important change in the flow dynamics.

 
The before mentioned observations are in agreement with the transient growth analysis performed by  Cantwell et al. <cit.> in the sense that the disturbance generated at the inlet grows while traveling downstream, mainly because of its interactions with the recirculation bubble. However, the authors would like to highlight some differences in our methodology and results when compared to <cit.>. Contrary to linear stability analysis with infinitesimal perturbations of <cit.>, DNS with finite amplitude perturbations are performed in the present work. In linear stability analysis, perturbations can either decay or grow, but in our case, perturbations could stay bounded for a very long time and then suddenly grow. The nonlinear interactions can also be observed in the DNS data, especially when there are interactions between the unsteady patterns and the recirculation bubble.










The WUP pattern can be also extracted from the flow by using a Reynolds averaging technique on the LS phase. The base flow is approximated by taking the average fields during the LS (from t=1500(s)  to t=1800(s)): 

    𝐮(x,y,z)=<𝐮(x,y,z,t)>_1500<t<1800.

Then the unsteady pattern can be extracted using:

    𝐮'(x,y,z,t)=𝐮(x,y,z,t)-𝐮(x,y,z).


Since the amplitude of 𝐮'(x,y,z,t) is small and stay bounded in laminar phase, it could be treated as a perturbation and qualitatively compared with linear stability analysis. To further compare our results with <cit.>, the coordinate system is converted  from Cartesian (𝐮'(x,y,z,t)) to the cylindrical (𝐮'(r,θ,z,t)). Then, the Fourier transform of 𝐮'(r,θ,z,t) in azimuthal direction, noted as 𝐮', is computed:

    𝐮'(r,k_θ,z,t)=1/2π∫_0^2π𝐮'(r,θ,z,t)e^-2iπθ k_θdθ,

where, the azimuthal modes energy are computed point-wise as:

    ( [ e_z(r,k_θ,z,t); e_y(r,k_θ,z,t); e_z(r,k_θ,z,t) ])=1/2( [ |u'_x(r,k_θ,z,t)|^2; |u'_y(r,k_θ,z,t)|^2; |u'_z(r,k_θ,z,t)|^2 ]),

and the total energy of each azimuthal modes is the sum of all three directions: 

    e_total(r,k_θ,z,t)=e_x(r,k_θ,z,t)+e_y(r,k_θ,z,t)+e_z(r,k_θ,z,t).









From the linear stability point of view, the perturbation 𝐮' could be decomposed into grow and decay modes. However, in our case, its magnitude stays bounded and fluctuates for a very long time, before the transient growth started and US2 emerge. It suggests there should be a very slow growing mode or new perturbation modes are introduced. This is what we investigated in this section. The observation of the time evolution of all azimuthal modes (over all z, r and k_θ) reveals  no steady increasing of energy. Instead, in addition to the regular WUP,  a perturbation appears in the steady zone (z<20d) and smoothly evolves into turbulence. 

A closer look into the time evolution of the perturbation amplitude in the vicinity of transition (1900<t<2200) is shown in figure <ref> for an arbitrary value of r (here r=0.25d) and three different streamwise positions. For clarity, only one value of r is shown, the other values of r exhibit similar behaviors. In this figure, the fluctuations of energy in streamwise direction (e_z(r,k_θ,z,t)) or spanwise direction (e_x(r,k_θ,z,t)+e_y(r,k_θ,z,t)) are recorded as a signal in time for given values of z, r and k_θ. From these signals, a local peak in the time evolution can be noticed that is located close to the transition moment. This peak can be recorded as a time instance t_p corresponding to the maximum of the energy signal. The value of collected t_p will be a function of spatial positions r, z and the mode k_θ. 

By comparing the fluctuation patterns, it is found that fluctuations in streamwise direction are around two orders of magnitude stronger than the sum of energy fluctuations in the two other spanwise directions in the laminar zone (first zone: z<20d). However, the gap in the energy levels gets closer at the transition point in the second zone: z≈50d and the non linear region, in the third zone: z>60d. 
Another observation is that the position of the peaks mainly depends on the streamwise position z and seems to be independent of spanwise variables (r and k_θ). For z≈ 17d, peaks start to have significant value  around t=2080(s). 

For a given position, the pattern of the peaks, u^p_z, can be extracted by subtracting the velocity fields at instance before the peak (t=2070(s)) from the instance at the peaks (t_p=2080 (s))

    u^p_z(x,y,z=17d)=u_z(x,y,z=17d,t=t_p=2080 (s))-u_z(x,y,z=17d,t=2070 (s)).



Figure <ref> presents the contour plots of different streamwise velocity patterns. It is found that the most intense perturbation in the peak pattern u^p_z (figure <ref>(a)) appears close to the strongest shear rate position in the streamwise velocity mean profile u_z (figure <ref>(b)). This suggests  the peak pattern is a shear instability and is independent of the regular unsteady pattern u'_z. The regular unsteady pattern (figure <ref>(c); see also equation <ref>) is the fluctuation collected before the peak emerge. The latter structure is completely different from both the peak pattern and the mean flow profile structures and the fact that it stays bounded for a large period of time suggests that it would not evolve into transition mechanism.











By fixing r at any arbitrary value (here we chose r=0.25d) and focusing on the first 4 modes (k_θ∈{0,1,2,3}), one could notice that the position of local peak evolves smoothly and linearly in the streamwise direction in the range of 10d<z<50d. This suggests that the peak is an perturbation pattern that appears in addition to the regular pattern, it gets carried downstream by the main flow and meanwhile amplified. In figure <ref>(a), the tracking of the local peaks in space (z) and time (t) are shown (all the value of r are superposed), where a linear fit of z/V_0+constant with V_0≈ 0.667 is found for the same region.

The error of the tracking process mainly comes from the the output frequency of the data and the fluctuation in z direction. All the peaks found in the range 10d<z<40d are contained within the linear fit with the error bar of ± 2δ_t where δ_t is the distance traveled by the peak between two consecutive outputs. It is interesting to mention that the peak evolution tracking in space (z) and in time (t) collapse for all the values of spanwise variables (r and k_θ). This indicates  the perturbation pattern is carried downstream by the same velocity V_0, for all ∀ r and ∀ k_θ. V_0 is expected to be close to the velocity profile where the perturbation is the most intense. The disturbance is generated punctually at (x_0,y_0,z_0), gets carried downstream with the mean velocity in the same position u_z(x_0,y_0,z_0) and then distributed in 2 spanwise directions. The maximum value of u^p_z(x,y) corresponds to x=-0.0097d and y=-0.1221d. At this position, the mean flow has the value of 0.58, which is roughly close to V_0.



When the pattern of u^p_z is carried out of zone 2 and enters zone 3 (at z≈50d), its velocity is slowed down. From this point, the non linearity starts to emerge, the peak is diffused and merge with another peaks that emerged earlier. As a consequence the localized peak in time is transformed into a plateau. In figure <ref>(b), the magnitude of the peaks over time are plotted for the first 4 azimuthal modes. The growth of all the modes follow a well defined power law, until z≈ 50 d where they start to saturate. 


Finally, we believe the chaotic behaviours that emerged in the third zone (z>50d) are the same convective instabilities that are generated by the shear in the first zone (z≈17d). The generation and evolution of the shear instability seems to be independent of the existence of other unsteady behaviours. It is noted that the phenomenon of the convective instability is also observed in the numerical simulation of flow with backward-facing step <cit.>.





§ HYSTERESIS















The hysteresis can appears when two flow solutions can exist for the same Re. Therefore, the initial conditions and mainly the parameters of the disturbance control the appearance of the two solutions. To test the hysteresis behaviour, two branches of the simulations with Re increasing and Re decreasing are investigated. The Re increasing branch starts from laminar flow, whereas the Re decreasing branch begins from an unsteady state, here is US2. This approach was also considered by  <cit.> using a tilt disturbance of amplitude δ. Depending on δ and Re a domain of hysteresis was observed by the same authors. Specifically, for δ=0.001, the coexistence region was reported for 1475<Re<1850. Moreover, they found that the hysteresis region grows as δ decreases. It is not possible to directly compare the effect of the tilt disturbance with our vortex disturbance because these perturbations are of different nature: the vortex disturbance perturbation introduces rotation, whereas the tilt disturbance introduces a translation to the flow that does not obey non-slip velocity condition at the wall and can therefore be considered unphysical. However, the following results discuss the universality of the hysteresis behaviour.

A series of eight simulations are performed, with a fixed amplitude of vortex perturbation, 𝒜=0.2. The initial condition is case 3a with Re= 1300. When the simulation time reached 1500s, the final state is analysed and used as the initial condition for the next run at a higher Re. Then, Re, is increased with step  of 25, up to Re=1400 (case 4a, 5a, 7a, 8a), as depicted in the space-time diagrams of figure <ref>(a). In figure <ref>(b), the decreasing Re branch is initialised with a laminar state (case 3a), then  Re is directly increased to 1400 (case 8b) and the decreasing path down to 1300 with a step of 25 (case 7b, 5b, 4b). The simulations from the increasing and decreasing paths are compared for the same Re the results are presented in figure <ref>(a) and <ref>(b) and show minor changes.
Looking at the drag coefficient also represented in figure <ref>(a) and (b), every change in Re initiates a peak or a transient increase of C_z corresponding to the turbulent spot that propagates and decays downstream. As seen no hysteresis is found while following this standard procedure.

Additional simulations using different Re steps and a larger range of Re, were performed. The decreasing branch started with Re=2000. Then Re, is decreased to 1700, 1350 and 1000 (simulation L1d, L2d, L3d). The results is compared with the increasing branch which start from Re=1000 then increase to 1350, then 1700 (simulation L3d, L2i, L1i). Again, the results are represented in the form of space-time diagrams, in figure <ref>(a) and <ref>(b). The data show minor differences in the space-time behaviour. However, this time the drag C_z, suggests a small loop of hysteresis. 

When comparing two series of simulation, one could notice: in the first loops, with small steps of Re (Δ Re=25), the two cases with Re=1350 are laminar in both increasing and decreasing branches. Whereas, in the second series, with larger steps of Re (Δ Re=350), the two cases with same Re=1350 show US2. The space-time diagram of 4 cases with Re=1350 are extracted in the two series and show in figure <ref>. The different behaviour of the same Re suggest that larger steps of Re could eventually trigger the unsteady behaviours and potentially hysteresis loop sooner than smaller steps. 


To quantify the extent of hysteresis behaviour, we first define the integral of C_z over Re for both increasing and decreasing branches, 

    S=∫ C_z(Re)dRe ,

 and then the relative difference is used to quantify the hysteresis behaviour:

    H=Δ S/S.

 
Here Δ S is the difference between 2 branches and S is the mean value. It is noted that in order to observe the hysteresis phenomenon clearly, a specific procedure needs to be implemented. In the decreasing branch, the variation of Re should be large enough to avoid the transformation from US1 to US2. On the other hand, in the increasing branch, the variation of Re should be small enough to keep  the flow laminar. Such loops is presented in figure <ref>, with the decreasing branch initiates at Re = 2000, which is decreased consecutively to 1700, 1350 and 1325 (simulation L1d, L2d, L3dbis). The increasing branch is initiated at Re =  1300 and is then increased consecutively to 1325, 1350 and 1700 (simulation 3a, 4a, 10a). Based on the criteria defined in equation <ref>, this procedure lead to a hysteresis of H=27.87% compare to the two previous loops with the measure of hysteresis are H=1.15% and H=0.75% respectively. 



In figure <ref>, the C_z from the last 100 seconds of each case is plotted against Re. The laminar states, Re= 1325 and Re = 1350, lead to C_z with almost the same value within 0.04%.  The unsteady state at Re = 1375 and Re = 1400 have slightly different final value of C_z within 3% because of the unsteady nature of the flow.
It is noted that the systematic study of all possible steps with extremely long time scales would be a tedious investigation thus beyond the scope of this study.




§ CONCLUSION

In the subcritical Re, range studied here, 1350≲ Re≲ 2000, it was shown that the flow through a sudden expansion in a circular pipe perturbed by a finite amplitude and a constant perturbation could be forced into unsteady states. With a strong enough inlet disturbance, that can break the recirculation region, the unsteady pattern gains the energy during the process and becomes sustainable even when the perturbation is removed. If the perturbation is not strong enough, the flow might return to its laminar state, and then develop an unsteady pattern with a scenario similar to  a linear growth process described by <cit.>, in which an infinitesimal upstream disturbance can be amplified while being carried downstream, when it reaches certain threshold, it can trigger the instability. This second type of instability exists only within a very narrow band of flow parameters, and it can start emerging at later time. However, it is a state that the flow must cross if it evolves gradually from a laminar state to a turbulence state. 


This marginal unsteady state has several distinguish feature compared to the fully developed turbulent state: it does not exhibit hysteresis behaviour at large observation time. With an appropriate parameters (small Re and 𝒜), it could decay and the flow can revisit the laminar state or it could self-sustain and the unsteady patch stay at the end of the recirculation zone, it generally has weaker amplitude than the first unsteady state. 
Either way, the time and space signal for either case, i.e. LS or US2, have some similarities. First they show a steady time-independent region right after the expansion at 0<z<20d; and second  they have a very well organised weakly unsteady patterns, so-called WUP, till the region before the reattachment zone (i.e. 20d<z<50d). However, after the reattachment zone, i.e z>50d, the low frequency disturbances can emerge and create a still bounded disturbance with larger wave length in the case of LS or can be amplified into a strongly time dependent pattern in the case of US2. The perturbation amplitude required for the first instability to emerge, seems to scale with the Reynolds number as a power law of -3, which is much steeper than the previous numerical estimation of <cit.>, who showed a power of -0.006.

Furthermore, it was found that the velocity fluctuations in the streamwise direction are dominant (i.e. about two order of magnitude larger) in the first two before-mentioned regions when compared to two other spanwise directions. However, this gap in the velocity fluctuation energy content gets closer after the transition point where the nonlinearities arise. Additionally, a peak fluctuation was observed which was independent of the spanwise variables. By comparing its most intense position to the shear rate map, it was shown that this peak pattern is a shear instability and appears on top of the regular unsteady pattern u'_z, where it carried downstream with constant speed and its amplitude grows exponentially (e^p_z ∼ z^n).

Finally, a hysteresis quantification procedure was established to clarify the ambiguities in the recent literature. It is found that for fairly small change, increasing or decreasing, in the Reynolds number and sufficiently large computational time, there exists no hysteresis unlike what is reported in <cit.>. However, if the change between the initial and the final quasi-steady flow is large enough, for instance by suddenly dropping the Reynolds number by a significant value, the turbulent flow inside the sudden expansion pipe can be self-sustained. For such cases, hysteresis behaviours are found within the reported physical time in the manuscript.



This work was granted access to HPC resources of IDRIS under the allocation 2017-100752 made by GENCI (Grand Equipement National de Calcul Intensif- A0022A10103).  The authors also acknowledge the access to HPC resources of French regional computing center of Normandy named CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie) under the allocations 2017002. The authors gratefully acknowledge financial support from LabEx “Laboratoires d’Excellence” project EMC3 “Energy Materials and Clean Combustion Center” through the project INTRA (LabEx EMC3-2016-PC ). Our work has also benefited from helpful discussions with Prof. A. P. Willis (University of Sheffield, UK), who suggested this particular form for the vortex perturbation.




jfm
