1. Introduction
In the context of continuous crystallisation processes, the authors have experimentally evidenced that introducing a flow restriction in a tubular micro-channel influences the nucleation rate, such as in the production of crystals of pharmaceutical ingredients (Rimez et al. Reference Rimez, Debuysschère, Conté, Lecomte-Norrant, Gourdon, Cognet and Scheid2018; Rimez, Debuysschère & Scheid Reference Rimez, Debuysschère and Scheid2019). The same kind of flow restriction has been used by Park, Song & Jung (Reference Park, Song and Jung2009) to sort particles in suspension, whereas Chang et al. (Reference Chang, Huang, Li and Jen2010) used a sudden expansion to sort blood components. Restrictions have also been used in heat exchangers (Zohir, Abdel Aziz & Habib Reference Zohir, Abdel Aziz and Habib2011), combustion systems and reactors (Hallett & Gunther Reference Hallett and Gunther1984) and in other microfluidic devices (Kadivar & Farrokhbin Reference Kadivar and Farrokhbin2017). Due to the appearance of such geometry in industries and laboratories, the behaviour of flow after a sudden expansion has already been widely studied as for pipe flow (3-D cylindrical) (Back & Roschke Reference Back and Roschke1972; Latornell & Pollard Reference Latornell and Pollard1986) than for channel flows (2-D and 3-D) (Tsai et al. Reference Tsai, Chen, Wang, Lin and Fu2006; Dagtekin & Unsal Reference Dagtekin and Unsal2011). It is well known that the shear layer after an expansion separates from the wall and creates a recirculation zone of length 
$X_{{r}}$ (figure 1a) before reattaching to the wall further downstream. For both pipe flow (Mullin et al. Reference Mullin, Seddon, Mantle and Sederman2009; Sanmiguel-Rojas & Mullin Reference Sanmiguel-Rojas and Mullin2012) and 2-D channel flow (Battaglia et al. Reference Battaglia, Tavener, Kulkarni and Merkle1997; Drikakis Reference Drikakis1997; Battaglia & Papadopoulos Reference Battaglia and Papadopoulos2005; Fani, Camarri & Salvetti Reference Fani, Camarri and Salvetti2012; Qian et al. Reference Qian, Ma, Zhan, Luo, Wang and Wang2016), it was found that increasing the Reynolds number (
${\textit {Re}}$) beyond a critical Reynolds value (
${\textit {Re}}_c$) induces a symmetry breaking of the recirculation zones developed after the expansion. This transition corresponds to a Pitchfork bifurcation. Realizing a linear stability analysis of such a system enables an efficient determination of the 
${\textit {Re}}_c$, and thus allows the effect of the expansion ratio (
$ER=D/d$), defined as the ratio between the downstream width (
$D$) and the upstream width (
$d$) of the sudden expansion, to be explored.
(a) Streamlines obtained for an expansion ratio: 
$ER={D}/{d}=3$ and 
${\textit {Re}}=100$ with the inlet velocity profile corresponding to a Poiseuille flow, 
$X_{{r}}$ corresponding to the recirculation length; (b) definition of the geometry, 
$d$ and 
$D$ the inlet and outlet widths of the channel, 
$L$ its length, 
$\ell$ the upstream length (§ 2.3). The grey area corresponds to the computational domain.
It has been shown that 
${\textit {Re}}_c$ is affected by the rheological behaviour of the fluid in the case of non-Newtonian fluids (Mishra & Jayaraman Reference Mishra and Jayaraman2002; Ternik Reference Ternik2009). In this study, we investigate another ingredient that might influence the stability of the flow after a 2-D sudden expansion, namely the shape of the inlet velocity profile. Indeed, it is known that for velocity profiles that deviate from a parabolic shape, an increase in the critical Reynolds number for the occurrence of turbulence in straight channel is observed (Potter Reference Potter1971). But to the best of our knowledge, the influence of the inlet velocity profile on critical Reynolds number for the symmetry breaking of the flow pattern in sudden expansions has never been studied. In this paper we focus on transition profiles, arising at the outlet of a local restriction (Durst & Loy Reference Durst and Loy1985), such as those used in our previous crystallization set-up (Rimez et al. Reference Rimez, Debuysschère, Conté, Lecomte-Norrant, Gourdon, Cognet and Scheid2018, Reference Rimez, Debuysschère and Scheid2019), even though only in a 2-D configuration.
The analyses realized in this study consist of considering various transition velocity profiles from plug to Poiseuille flows and simulating with direct numerical simulation (DNS) the behaviour of the recirculation length 
$X_{{r}}$ for several 
$ER$ and 
${\textit {Re}}$ in order to determine the 
${\textit {Re}}_c$ for each condition. Linear stability analysis (LSA) has been realized using various inlet profiles to validate the results obtained with time dependent simulation and to extend the parametrical analysis at lower computational cost. Finally, the length of the recirculation zone has been correlated to the convection–diffusion transport of momentum.
2. Methodology
The geometry of the problem along with the Cartesian coordinate system 
$(x',y')$ are given in figure 1. For the two approaches, DNS and LSA, the variables are made dimensionless using 
$u_{{mean}}$, the mean velocity, 
$\rho u^2_{{mean}}$, the pressure scale where 
$\rho$ is the density, 
$d$ the upstream width of the 2-D channel and 
$d/u_{{mean}}$ the timescale. The meshes used for the modelling are available in § 2.4, and convergence tests have been performed to validate them.
2.1. Direct numerical simulation
The 2-D fluid behaviour is modelled using the dimensionless Navier–Stokes and continuity equations for incompressible fluids:
\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u}+\boldsymbol{\nabla} p-\frac{1}{{\textit{Re}}}{\rm \Delta}\boldsymbol{u}=0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0, \end{gather}
where 
$\boldsymbol {u}=(u,v)$ with 
$u$ and 
$v$ the velocities in the 
$x$ and 
$y$ directions, 
$t$ the time, 
$p$ the pressure field and 
${\textit {Re}} = u_{{mean}}d/\nu$ the Reynolds number with 
$\nu$ the kinematic viscosity. A no-slip boundary condition is imposed at the walls
\begin{equation} \boldsymbol{u}=0, \end{equation}
the inlet velocity is set as function of 
$y$,
\begin{equation} u(0,y)= f(y) \quad\forall y \in [-1/2 < y < 1/2] \end{equation}
and a free outflow is set at the outlet, 
$x=\bar {L}$, with 
$\bar {L}=L/d$,
\begin{equation} p-\frac{1}{{\textit{Re}}}\partial_x u=0\quad \mbox{and} \quad \partial_x v=0. \end{equation}
These equations are solved using the open source DNS code Nek5000 based on the spectral elements method (SEM) and using seventh-order polynomials (P7–P5) discretization (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). Because we are looking for stationary solutions, the convergence criterion to stop the simulations is 
${||u_{t}-u_{t-{\rm \Delta} t}||}/{{\rm \Delta} t}<10^{-7}$, where 
${\rm \Delta} t$ is the time step. The dimensionless length of a recirculation zone, 
$\bar {X}_{{r}}$, is defined as the downstream distance from 
$x=0$ at which the wall shear stress 
$\partial _y u$ is zero. In order to capture the asymmetry, for each expansion ratio, the procedure is two-fold: (i) obtain a stationary asymmetric solution by choosing a value of 
${\textit {Re}}$ sufficiently large to be beyond 
${\textit {Re}}_c$ and by imposing as initial condition a forcing of sufficiently large amplitude to ensure the asymmetry; (ii) use the preceding solution as an initial solution and run simulations by decreasing the Reynolds number until symmetric solutions are recovered, indicating that 
${\textit {Re}}_c$ has been reached. 
${\textit {Re}}_c$ has been defined as the average Reynolds number within the interval for which a symmetry breaking was observed.
2.2. Linear stability analysis
Following classical linear stability theory, the flow variables (
$\boldsymbol {u},p$) are characterized by their steady base flow conterparts (
$\boldsymbol {u}_{{b}},p_b$) such that
\begin{gather} (\boldsymbol{u}_{{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u}_{{b}}+\boldsymbol{\nabla} p_{{b}}-\frac{1}{{\textit{Re}}}{\rm \Delta}\boldsymbol{u}_{{b}}=0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_{{b}} =0, \end{gather}closed by the same boundary conditions (2.3)–(2.5a,b). The flow is also characterized by its unsteady perturbed part and can therefore be decomposed as
\begin{gather} \boldsymbol{u}(x,y,t)=\boldsymbol{u}_{{b}}(x,y)+\epsilon \tilde{\boldsymbol{u}}(x,y)\,\textrm{e}^{\sigma t}, \end{gather}
\begin{gather}p(x,y,t)=p_{{b}}(x,y)+\epsilon \tilde{p}(x,y)\,\textrm{e}^{\sigma t}, \end{gather}
where 
$\epsilon \ll 1$ is a small order parameter, 
$\tilde {\boldsymbol {u}}$ and 
$\tilde {p}$ are the complex perturbations and 
$\sigma$ is the complex eigenvalue. By substituting (2.8) and (2.9) into (2.1) and (2.2) and linearizing with respect to 
$\epsilon$, we obtain
\begin{gather} \sigma \tilde{\boldsymbol{u}}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u}_{{b}}+ (\boldsymbol{u}_{{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \tilde{\boldsymbol{u}} +\boldsymbol{\nabla}\tilde{p} -\frac{1}{{\textit{Re}}}{\rm \Delta}\tilde{\boldsymbol{u}}=0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} =0, \end{gather}together with the inlet and walls boundary conditions,
\begin{equation} \tilde{\boldsymbol{u}}=0 , \end{equation}and the outlet condition,
\begin{equation} \tilde{p}-\frac{1}{{\textit{Re}}}\partial_x \tilde{u}=0\quad\mbox{and}\quad\partial_x \tilde{v}=0. \end{equation}
The problem is an eigenvalue problem where 
$\sigma =\lambda +\textrm {i}\omega$, where 
$\lambda$ is the growth rate and 
$\omega$ is the pulsation. Equations are discretized by a finite element method (FEM) using P2–P1 Taylor–Hood elements. Meshes and matrices are generated with the software FreeFem++. The base-flow described by (2.6)–(2.7) is solved with no symmetry condition with a Newton–Raphson iterative method, and the solution of the linear system is computed via the Unsymmetric Multifrontal sparse LU Factorization PACKage (UMFPACK) (Fani et al. Reference Fani, Camarri and Salvetti2012; Hecht Reference Hecht2012). The flow is linearly stable when the real parts of all eigenvalues are negative (
$\lambda <0$) and unstable when at least one of them is positive (
$\lambda >0$). The leading global mode, defined as the global mode with the largest growth rate, is investigated to determine the global stability of the base flow. Figure 2(a) shows the growth rate increases of the leading global mode while the Reynolds number increases for an 
$ER$ of 
$15$ with a plug flow as inlet velocity. Similar to § 2.1, the critical Reynolds number (
${\textit {Re}}_c$) is, here, defined as the average Reynolds number in the interval where the real part of the leading eigenvalue is positive. For example, in figure 2(a), the growth rate of the leading mode for 
${\textit {Re}}=16.050$ is negative and is positive for 
${\textit {Re}}=16.083$; the critical Reynolds number is then set as 
${\textit {Re}}_c=16.067\pm 0.017$. Figure 2(b) depicts the streamwise (
$\tilde {u}$) and the cross-stream (
$\tilde {v}$) velocity perturbations and the pressure perturbation (
$\tilde {p}$) of the leading global mode at 
${\textit {Re}}=20$. The streamwise velocity perturbation and the pressure perturbation are anti-symmetric regarding the 
$y$ axis, and the cross-stream velocity is symmetric. In all the following results (see § 3), we noticed neither the appearance of a second unstable global mode nor the appearance of non-zero pulsation.
(a) Growth rate (
$\lambda$) of the leading global mode obtained with the LSA for 
$ER=15$ with a plug flow as inlet velocity; (b) perturbation for 
${\textit {Re}}=20$ for the velocity in the streamwise direction (
$\tilde {u}$), in the cross-stream direction (
$\tilde {v}$) and the pressure field (
$\tilde {p}$).
2.3. Parametrization of the inlet velocity profile
Since we are interested in studying the influence of the inlet velocity profile on the flow stability in a 2-D sudden expansion, we first model the flow in a straight channel of finite length 
$\ell$ where a plug flow is imposed at the inlet. Thereafter, we equalized the convection time with the diffusion time:
\begin{equation} \frac{\ell}{u_{{mean}}}\simeq\frac{d^2}{\nu}\rightarrow \frac{\ell}{d} \simeq \frac{u_{{mean}}d}{\nu} = {\textit{Re}}. \end{equation}
From this equality, the dimensionless channel length 
$\ell /d$, corresponding to the development length of the velocity profile, can be rescaled with the Reynolds number and defined as 
$k=\ell /(d{\textit {Re}})$. This parameter means that for two different Reynolds numbers (
${\textit {Re}}_{1}$ and 
${\textit {Re}}_2$), the flow profiles are equivalent while 
$\ell _{1}/(\textrm {d} {\textit {Re}}_1)=\ell _2/(d {\textit {Re}}_2)$. The validity of this hypothesis (see appendix A) is shown to be satisfied for 
${\textit {Re}} \gtrsim 333$. As a result, for the following analysis, the profiles obtained at 
${\textit {Re}}=400$ are represented in figure 3. These profiles are then imposed at the inlet of the sudden expansion such as schematized in figure 1. Even though the theoretical 
$k$ value for a Poiseuille flow should correspond to infinity, we used the value 
$k=0.045$, which corresponds to the rescaled dimensionless distance from the inlet at which the discrepancy between the central velocity and the theoretical central velocity of a Poiseuille velocity profile is less than 1 % to parametrize (3.2). This 
$k$ value is comparable to the one obtained in the literature (Lautrup Reference Lautrup2011).
Velocity profiles at the rescaled downstream development length 
$k=\ell /(d{\textit {Re}})$ in a straight channel for 
${\textit {Re}}=400$.
2.4. Meshes
For both DNS and LSA, the computational domains are built with an inlet of width 1, an outlet of width 
$ER$ and a length of 
$50$ – i.e. long enough to ensure that the results are independent of this length.
The DNS mesh has 
$120$ elements along the flow direction with a ratio of 
$1.01$. In the cross-flow direction, the channel is divided in three parts where the meshes are built according to the Chebyshev nodes theory ((2.15), where 
$y_i$ is the position of the 
$i{{\rm th}}$ node). The central zone, which has the same width as the inlet, contains 11 elements. The two other zones at each side of the central one each contain 16 elements and have a width 
$(ER-1)/2$.
\begin{equation} y_i=\frac{1}{2}(a+b)+\frac{1}{2}(b-a)\cos\left(\frac{2i-1}{2n}\right);\quad (i=1,\ldots,n). \end{equation}
The mesh used for the LSA is described in table 1, where the positions (
$x_i,y_i$) of the 
$i{{\rm th}}$ node are described as a function of 
$ER$ and a constant 
$a$ varying with 
$ER$. For the two analyses, DNS and LSA, the convergence of the meshes have been validated by increasing the number of elements which induced no differences, respectively, for the length of the recirculation zones and for the eignevalues.
Left: Mesh used for the LSA, position (
$x_i,y_i$) of the 
$i{\rm th}$ node for 
$i$ between 
$0$ and 
$n$ the total number of node in one direction. Right: value of the constant 
$a$ for the different expansion ratios 
$ER$.
3. Results
3.1. Flow behaviour
The DNS modelling was carried out for different profiles characterized by their 
$k$ parameters, including the two limit cases (Poiseuille and plug flows) and for two values of 
$ER=2$ and 
$ER=3$. Figure 4 shows Pitchfork bifurcations, indicating that the dimensionless recirculation length (
$\bar {X}_r$) at both sides of the channel do not have the same size when the Reynolds number is above the critical value, 
${\textit {Re}}_c$. The two cases (a) and (b) show an increase of 
${\textit {Re}}_c$ with the decrease of 
$k$. It can be seen that the slope before the supercritical Pitchfork bifurcation decreases when 
$k$ is decreased, meaning that at a fixed Reynolds number, the reattachment length 
$\bar {X}_{{r}}$ is always shorter for a plug flow than for a Poiseuille flow at the inlet.
Pitchfork bifurcations obtained by DNS for the length (
$\bar {X}_{{r}}$) of the recirculation zones for (a) 
$ER=2$ and (b) 
$ER=3$ for inlet velocity profiles going from plug flow to Poiseuille flow, as parametrized by 
$k$. Dashed lines are the unstable recirculation lengths for the two extreme profiles.
3.2. Critical Reynolds number
On extending the LSA to other 
$ER$ (2.4, 2.6, 3.5, 4, 6, 8, 15) and profiles (
$k=0.00632$ and 
$k=0.03304$), figure 5(a) shows that the critical Reynolds numbers decrease when 
$ER$ is increased (see appendix B for the raw data). Comparing the critical Reynolds numbers obtained in the linear stability analysis at 
$ER=2$ and 
$ER=3$ with the ones obtained with DNS (figure 4) indicates that the results are in agreement. To rationalize these results, we define a 
$\delta (k)$ factor in order to correlate the 
${\textit {Re}}_c(ER,k)$ for a given 
$ER$ and for a known 
$k$ profile with the critical Reynolds numbers for a plug flow (
${\textit {Re}}_c^{Pl}$) and a Poiseuille flow (
${\textit {Re}}_c^{Po}$) at the inlet:
\begin{equation} {\textit{Re}}_c(ER,k)\simeq\delta(k)[{\textit{Re}}_c^{Pl}(ER)-{\textit{Re}}_c^{Po}(ER)]+{\textit{Re}}^{Po}_c(ER). \end{equation}(a) Critical Reynolds numbers as a function of the expansion ratio 
$ER$, calculated with the LSA for inlet velocity profiles varying from a plug flow to a Poiseuille flow; (b) values of the 
$\delta$ factor as function of the rescaled dimensionless distance 
$k$. The dotted line is a fitting by the least square method (see text).
By plotting the function 
$\delta (k)$ in figure 5(b), a smooth decrease of the value of 
$\delta$ from 
$1$ (plug flow) to 
$0$ (Poiseuille flow) is observed for increasing 
$k$. This relation has been fitted using the least squares method and by imposing two cross points at 
$k=0$ (
$\delta =1$) and 
$k=0.045$ (
$\delta =0$) corresponding to the two limit profiles, which gives
\begin{equation} \delta = \frac{\beta}{k^n+\gamma}+\zeta, \end{equation}
where 
$n=0.785$, 
$\gamma =0.012$, 
$\beta =(\gamma\, 0.03^n+\gamma ^2)/(0.03^n)= 0.014$ and 
$\zeta =1-\beta /\gamma =-0.136$. Such a curve is useful to predict the critical Reynolds number for each rescaled dimensionless 
$k$ and 
$ER$, by using only the 
${\textit {Re}}_c$ for the two limit profiles (Poiseuille and plug flows), which are given in table 2.
Critical Reynolds numbers for the Pitchfork bifurcation depending on the inlet velocity profile (Poiseuille or plug flows) and the expansion ratio 
$ER$.
It is interesting to note that the critical Reynolds numbers can also be approximated as a function of the expansion ratio and the inlet flow profile,
\begin{equation} {\textit{Re}}_c=\kappa \frac{1+\delta}{ER^m}+\xi, \end{equation}
where 
$m=3.06$, 
$\kappa =1298.39$ and 
$\xi =15.74$ are fitting parameters. This correlation has been drawn in figure 5(a) for the two inlet limit cases, namely Poiseuille and plug flow profiles, for which 
$\delta$ is respectively equal to zero and one. Interestingly, when one extrapolates the values for the plug flow (
$\delta =1$) at 
$ER=1$, it appears that the value of 
${\textit {Re}}_c = 2613$ is close to 
$2700$, i.e. the value for the laminar to turbulent transition in a straight channel as referenced in the literature (Morini Reference Morini2004; Wang et al. Reference Wang, Gao, Tan and Xu2012).
3.3. Scaling for the length of the recirculation zones
Figure 4 shows that the slope of the recirculation length below the critical Reynolds numbers depends on the profile. While extending the computation to 
$ER=4,6$ and 
$8$ by DNS, the same trends have been obtained (see figure 6a). With the aim of finding a scaling law for this observation, we then proceed by assuming that the slope variations can be rationalized using a convection–diffusion argument. For that purpose, we introduce the non-dimensional parameter 
$\alpha$, a function of 
$k$, which describes a fraction of the distance 
$d$ in addition to the step height (
$h$) through which lateral diffusion of momentum has to take place as schematized in figure 7 in dimensionless form. The associated time scale can thus be written as
\begin{equation} t_\nu \simeq\frac{(h+\alpha d)^2}{\nu}. \end{equation}
By inspection of the velocity profiles, we have found that 
$\alpha$ is correlated to the previously set parameter 
$\delta$ by 
$\delta =1-2\alpha$. Therefore, 
$\delta$ corresponds to the central non-dimensional distance (see figure 7) over which the velocity profile is almost flat (see figure 8b), while 
$\alpha$ characterizes the complementary non-dimensional distance over which the velocity variations take place. Equalizing this lateral diffusion time scale and the convection time regarding the development length (
$X_{{r}}$) of the recirculation, yields
\begin{equation} \frac{X_{{r}}}{u_{{mean}}}\propto\frac{\left[h+\dfrac{(1-\delta)}{2} d\right]^2}{\nu}, \end{equation}
or in dimensionless form and by using 
$h=(D-d)/2$,
\begin{equation} \bar{X}_{{r}}\propto {\textit{Re}} \left(ER-\delta \right)^2={\textit{Re}}^{*}. \end{equation}
Using the 
$\delta$ factors established in § 3.2, both the recirculation lengths obtained before the Pitchfork bifurcations (figure 6b) and their slopes (inset of figure 6a) are seen to collapse for all the 
$ER$ and entry profiles studied. It must be noted that this collapse of data points in figure 6(b) does not go through 
$(0,0)$, which is in good agreement with the results of Khodaparast, Borhani & Thome (Reference Khodaparast, Borhani and Thome2014), which show the nonlinear behaviour of the recirculation length at low Reynolds numbers, due to the effect of the non-negligible streamwise diffusion of momentum (see also appendix A).
(a) Recirculation length obtained before the Pitchfork bifurcation for 
$ER = 2$, 3, 4, 6 and 8 for different inlet velocity profiles going from plug to Poiseuille flow profiles; the inset depicts the square root of the slopes as a function of the expansion ratio and the 
$\delta$ parameter. (b) Collapse of the recirculation lengths as a function of the rescaled Reynolds, 
$Re^*$ obtained in (3.6).
Definition of the dimensionless parameters 
$\alpha$ and 
$\delta$ influencing the establishment of the recirculation zone where 
$\bar {X}_r$ and 
$\bar {h}$ are the dimensionless recirculation length and height, and 
$ER$ is the expansion ratio.
(a) The 
$y$-derivative of the streamwise velocity component (b) in black-dashed region of size 
$\delta$ (schematized in red) where the momentum diffusion is considered as negligible.
The physical interpretation of 
$\delta$ as the extension of the shear-free central jet region is corroborated by plotting the position 
$y=\delta /2$ on each derivative of the studied inlet velocity profile with a black point in figure 8(a). These black points all fall within a small grey zone corresponding to 
$|{\partial {u}}/{\partial {y}}|<0.5$, meaning that the shear in the interval 
$[{-\delta }/{2},{\delta }/{2}]$ is close to zero. Figure 8(b) represents these inlet velocity profiles where the black-dashed part of dimensionless length 
$\delta (k)$ represents the zone where momentum diffusion is assumed to be negligible. At both sides of the shear-free central zone, the shear is large, and the associated vorticity is assumed to actively participate in the formation of the recirculation zones.
4. Discussion and conclusions
Flow recirculations were studied downstream of a sudden expansion in a 2-D channel. Those flow behaviours were analysed by DNS and by means of LSA. The two techniques were applied on a set of inlet velocity profiles ranging from Poiseuille flow to plug flow as limits. The intermediate profiles were parametrized by a 
$k$ factor, representing the rescaled dimensionless distance downstream of the inlet in a straight channel where a plug flow was imposed at the inlet. Parametrized inlet profiles were then applied to different expansion ratios ranging from 
$2$ to 
$15$.
First, we studied the appearance of Pitchfork bifurcations. We noticed for a decreasing 
$k$ value an increase of the critical Reynolds number at which the Pitchfork bifurcation occurs. A correlation factor, 
$\delta$, has been introduced in order to predict 
${\textit {Re}}_c$ as a function of the 
$k$ constant and the expansion ratio. We then showed that the critical Reynolds numbers can be approximated by means of a function of 
$\delta$ and 
$ER$.
Second, we pointed out that the recirculation lengths were influenced by both the inlet flow profiles and by the expansion ratios. The length of the recirculation zones tends to increase as the flow profile gets closer to the Poiseuille flow profile. This observation has been correlated to momentum diffusion by means of the previously determined dimensionless factor 
$\delta (k)$. This separates the inlet profile into two zones: (i) the shear-free central region of extension 
$\delta (k)$, which has a negligible effect on the recirculation zones; (ii) the two regions at either side of the central one where the shear and the associated vorticity are significant, which then influence the recirculation zones. The introduction of this shear-free region 
$\delta (k)$ into the convection/diffusion equilibrium allowed us to point out a linear correlation between the length of the recirculation zones and the rescaled Reynolds number 
${\textit {Re}}^*$, taking the momentum diffusion into account.
Having demonstrated the strong influence of the inlet velocity profile on the flow structure and stability in a 2-D sudden expansion, it could be relevant to extend these analyses for the design of 3-D microdevices. However, while the length of the recirculation zone obtained in the present paper for any ER and constriction length is expected to be applicable for a 2-D axisymmetric flow, the critical Reynolds number for the Pitchfork bifuraction is not. Sanmiguel-Rojas & Mullin (Reference Sanmiguel-Rojas and Mullin2012) and Mullin et al. (Reference Mullin, Seddon, Mantle and Sederman2009) demonstrated that 
${\textit {Re}}_c$ for a pipe flow with a Poiseuille velocity profile at the inlet is much larger (
$ER=2$, 
${\textit {Re}}_c\sim 1139$) than the one for a 2-D configuration. They also showed that the critical Reynolds number is subject to a hysteresis phenomenon. Yet it is believed that having a non-developed flow at the inlet of a pipe with a sudden expansion should still have a stabilizing influence for the same reasons as those presented here for a 2-D geometry. In 3-D geometry, other parameters may also influence the flow behaviour. Vaniershot & Van den Bulck (Reference Vaniershot and Van den Bulck2008) demonstrated that the recirculation lengths are also influenced by upstream swirl for high Reynolds numbers, which could possibly be correlated to the modification of the inlet velocity conditions. At the same time, Cantwell, Barkley & Blackburn (Reference Cantwell, Barkley and Blackburn2010) demonstrated that sudden expansions cause the amplification of the inlet noise and Zhang & Luo (Reference Zhang and Luo2018) demonstrated that sudden constrictions may induce an increase of the swirl intensity, which therefore promotes mixing. All these parameters should be in the scope of future studies and may eventually lead to help the design of the following: (i) sorting microsystems such as those described by Kadivar & Farrokhbin (Reference Kadivar and Farrokhbin2017), Park et al. (Reference Park, Song and Jung2009), Chang et al. (Reference Chang, Huang, Li and Jen2010) or Volpe, Gaudiuso & Ancona (Reference Volpe, Gaudiuso and Ancona2019), which may possibly be improved by controlling the size of the recirculation zones; (ii) microreactors such as the microcrystallizer described by (Rimez et al. Reference Rimez, Debuysschère and Scheid2019), which may also be improved by promoting passive mixing.
Acknowledgements
This work was executed thanks to the Walloon Region and its financial support during the MecaTech-Legomedic project. R.D. and B.S. thank the F.R.S.-FNRS for financial support.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Profile convergence
Based on the simulations described in § 2.3, velocity profiles have been extracted for Reynolds numbers from 33.3 to 533.3 at different rescaled downstream distances (
$k$) from the inlet. Figures 9(a)–9(d) show that profiles are not equivalent for an increasing Reynolds number. One observes that the profiles collapse with 
$Re\gtrsim 333.3$. Thereby, we arbitrarily choose for reference the profiles obtained with a Reynolds number of 
$400$. Thanks to the least squares method, for various Reynolds numbers, we compare in figure 10 the downstream distance k from the inlet at which the velocity profile is extracted with the one (
$k|_{Re=400}$) in the case of the chosen reference (
$Re=400$). The curves collapse when 
$Re\gtrsim 333.3$, confirming the first observation made in figure 9. This discrepancy within the downstream distances for low Reynolds numbers originate from the longitudinal diffusion of momentum. Indeed, equalizing the convection time and the longitudinal diffusion time gives
\begin{equation} \frac{\partial{\ell}}{{u_{mean}}}\simeq\frac{\partial{\ell}^2}{\nu}, \end{equation}
where 
$\partial {\ell }$ is the longitudinal diffusion length. It allows the dimensionless diffusion length (
$\partial {\ell }^*$) to be described as a function of the Reynolds number
\begin{equation} \partial{\ell}^*=\frac{\partial{\ell}}{d}\simeq\frac{\nu}{d{u_{mean}}}=\frac{1}{Re}. \end{equation}
For a given Reynolds number and for a velocity profile obtained at the dimensionless distance 
$k$ from the inlet, the dimensionless distance 
$\partial {\ell }^*$ can be approximated by the distance to recover the same velocity profile in the case of the reference Reynolds number (
$Re=400$),
\begin{equation} \partial{\ell}^*\propto k-k|_{Re=400}. \end{equation}
Figure 11 shows this length as a function of the Reynolds number for a profile extracted at 
$k|_{Re=400}=0.01522$. It appears that this is linearly correlated to 
$1/Re$ as predicted by (A 2), which confirms that the observed changes were due to longitudinal diffusion.
Computed velocity profiles by DNS at 
$k=0.00038$ (a), 
$k=0.0078$ (b) 
$k=0.01151$ (c) and 
$k=0.01522$ (d) for 
$Re$ between 33.3 and 533.3 in straight channel.
For various Reynolds numbers ranging from 33.3 to 533.3, downstream distances 
$k$ from the inlet at which the same profile is found as at a distance 
$k|_{Re=400}$ from the reference Reynolds number (
$Re=400$).
Correlation between 
$\partial {\ell }^*$ and Reynolds numbers for a velocity profile obtained at the downstream dimensionless distance from the inlet 
$k|_{Re=400}=0.01522$.
Appendix B. Data
Table 3 contains the numerical values of the critical Reynolds numbers (
$Re_c$) for the analysed expansion ratios (
$ER$) for various inlet velocity profiles ranging from plug flow to Poiseuille flow as parametrized by 
$k$. The results shown have been obtained by DNS or LSA. The values 
$\pm \varDelta$ describe the half step size used to determine 
$Re_c$.
Critical Reynolds numbers (
${\textit {Re}}_c$) obtained with DNS and LSA, and the half step sizes used to determine these values.
