1. Introduction
Sudden expansion geometries are used in several engineering and biomedical applications, such as augmenting heat transfer in heat exchangers (Zohir, Aziz & Habib Reference Zohir, Aziz and Habib2011), combustion chambers (Reddy, Sujith & Chakravarthy Reference Reddy, Sujith and Chakravarthy2006), fuel injectors, sorting blood cells in microfluidic devices (Chang et al. Reference Chang, Huang, Li and Jen2010), etc. It has also been used for increasing the nucleation rate of pharmaceuticals by an order of magnitude (Rimez, Debuysschère & Scheid Reference Rimez, Debuysschère and Scheid2020) and fostering the mixing rate in fluids (Lü et al. Reference Lü, Zhu, Wang and Luo2016). Nature's refined designs of expansion are present in the aortic root of the heart, arterial bifurcation, etc.
Early experimental works quantified the velocity patterns downstream of the sudden expansion (Durst, Melling & Whitelaw Reference Durst, Melling and Whitelaw1974; Cherdron, Durst & Whitelaw Reference Cherdron, Durst and Whitelaw1978; Fearn, Mullin & Cliffe Reference Fearn, Mullin and Cliffe1990; Durst, Pereira & Tropea Reference Durst, Pereira and Tropea1993). They reported a pair of recirculating eddies at the top and bottom walls after the expansion. At sufficiently low Reynolds numbers (Re), the recirculating eddies are equal in size. As Re is increased, this symmetry vanishes, with one of the recirculation regions growing in size at the expense of the other. Cherdron et al. (Reference Cherdron, Durst and Whitelaw1978) ascertained that the amplification of disturbances generated at the lee of the expansion by the shear layer is responsible for the asymmetric nature of flow above a critical Reynolds number. Fearn et al. (Reference Fearn, Mullin and Cliffe1990) and Durst et al. (Reference Durst, Pereira and Tropea1993) quantified the symmetry breaking of the problem numerically by solving a steady, two-dimensional Navier–Stokes equation and reported a bifurcation diagram. They obtained a critical Reynolds number, 
${{Re_{c}}}$, below which the symmetric solution is stable and above which two solutions that are antisymmetric with respect to each other exist. Sobey & Drazin (Reference Sobey and Drazin1986) categorized this instability as pitchfork bifurcation. The observed asymmetry is attributed to Coanda effects (Pitton, Quaini & Rozza Reference Pitton, Quaini and Rozza2017). Investigators have performed bifurcation analysis, linear stability analysis (LSA) and weakly nonlinear studies to determine the onset of primary and higher-order instabilities (Alleborn et al. Reference Alleborn, Nandakumar, Raszillier and Durst1997; Hawa & Rusak Reference Hawa and Rusak2001). Hawa & Rusak (Reference Hawa and Rusak2001) ascertained that the symmetric state can lose stability due to the combined effects of convection of base flow vorticity and its perturbation, respectively, by the perturbed and base axial flow along with viscous dissipation of vorticity perturbation. The authors also report an additional recirculation region at the side of the smaller eddy as Re increases.
Battaglia et al. (Reference Battaglia, Tavener, Kulkarni and Merkle1997) and Drikakis (Reference Drikakis1997), amongst others, have studied the effect of expansion ratio 
${{ER}}$, defined as 
${{ER}}=d/D$, where 
$d$ and 
$D$ are the heights of the inlet and outlet channel, respectively, on the bifurcation phenomena in the expansion channel. They found that as 
${{ER}}$ increases, 
${{Re_{c}}}$ decreases. The effect of the angle of expansion has also been studied by investigators (Shapira, Degani & Weihs Reference Shapira, Degani and Weihs1990; Jotkar & Govindarajan Reference Jotkar and Govindarajan2019). It was found that as the angle of expansion increased, the size of the recirculation eddy also increased, with 
${{Re_{c}}}$ dropping significantly. Drikakis (Reference Drikakis1997) also examined various numerical schemes in their ability to predict 
${{Re_{c}}}$ and found that higher-order finite-difference schemes accurately predict symmetry-breaking phenomena. Recently, Debuysschère et al. (Reference Debuysschère, Siconolfi, Rimez, Gallaire and Scheid2021) investigated the influence of inlet velocity profiles on flow stability in a sudden expansion channel. Their inlet profiles varied from plug to Poiseuille profiles, with several intermediate ones. They obtained higher 
${{Re_{c}}}$ with the plug than Poiseuille flow at the inlet, with an inverse relationship with the recirculation eddy length 
${{X_{r}}}$.
Most investigations on the stability of sudden expansion flow are restricted to steady inflow cases. In the realm of pulsatile inflow in sudden expansion, Ma, Li & Ku (Reference Ma, Li and Ku1994) have compared the role of steady and pulsatile inflow in heat transfer augmentation. They considered only a single time period of 8.7 s for the case of pulsating inflow and focused mainly on the heat and mass transfer aspects at different Prandtl numbers. Stephanoff et al. (Reference Stephanoff, Pedley, Lawrence and Secomb1983) considered a pulsating inflow in a one-sided smooth-expansion channel via altering the height of the channel. They obtained different flow regimes for low and moderate Strouhal numbers, 
$St=fd/\nu$, depending on the frequency of oscillation, channel height and mean velocity at the inlet. At lower 
$St$, the flow behaved quasi-steadily with a single recirculation eddy behind the expansion. At moderate 
$St$, however, a vortex wave-like phenomenon occurred with the presence of many small asymmetrically placed eddies at both walls and the core flow moving in a sinuous manner. Sobey (Reference Sobey1985) also observed this phenomenon in their expansion channel and attributed the cause of this sinuous motion to a shear-layer instability rather than a vortex wave.
Pulsatile sudden expansion flows are often encountered in arteries with stenosis. Neofytou & Drikakis (Reference Neofytou and Drikakis2003a) performed simulations for a two-dimensional, pulsatile flow of Newtonian and non-Newtonian fluids (Casson, power law and Quemada models) over one-sided stenotic blockage in a channel. Neofytou & Drikakis (Reference Neofytou and Drikakis2003b) observed that in steady flows over sudden expansion, both Newtonian and non-Newtonian fluids disembark symmetry-breaking bifurcations. Wu et al. (Reference Wu, Aubry, Antaki and Massoudi2020) reported that the recirculation zone after a sudden expansion arterial model with pulsatile inflow shows a reduced red blood cell concentration. They suggested that these recirculation zones may favour platelet accumulation. They found a red blood cell depletion region just after the sudden expansion, which might favour platelet accumulation. For high Reynolds number flows, Devenport & Sutton (Reference Devenport and Sutton1993) performed an experimental study of flow through axisymmetric sudden expansion with an expansion ratio of 1.875 at 
$Re=35\,000$. They found a turbulent free shear layer (similar to that of a round jet) due to the flow separation glances at the wall at an angle that reduces in the presence of a centrebody. In the latter case, higher curvature of the shear layer results in suppression of the turbulence levels. Karantonis et al. (Reference Karantonis, Kokkinakis, Thornber and Drikakis2021) studied compressibility effects on flows in a sudden expansion channel at high Reynolds number and moderate Mach numbers. They found that compressibility effects substantially influence the flow.
However, it may be observed that the effect of inflow pulsatility on the flow structures in sudden expansion channels has been scarcely reported, although it has relevance in several biological and engineering applications. Moreover, to the best of our knowledge, the symmetry-breaking characteristics due to a pulsatile inflow over an expansion have never been studied. Hence, in the present work we investigate the effect on flow stability through planar sudden expansion via inflow pulsatility. We vary the amplitude and period of pulsation to unveil different flow features. We show streamline patterns to qualitatively describe the flow and define parameters such as scaled recirculation lengths, 
${{X_{r}^{*}}}$, and scaled position of the centre of the recirculation region, 
${{X_{c}^{*}}}$, to quantify the size and location of the eddying structures. We show four kinds of flow patterns or modes depending on pulsation amplitude and frequency, which influence the 
${{Re_{c}}}$. Here, we denote the temporal evolution of recirculation zones as modes. We also report the influence of 
${{ER}}$ on stability.
Arterial flow is pulsatile in nature, and sudden increases in lumen area are often observed in the vascular networks (e.g. aortic or carotid sinus). Stenotic blockages also form in the arterial lumen due to cholesterol deposition. The stenosis portion resembles a smooth contraction–expansion. Depending on the site of the stenosis, the Reynolds number and Womersley number parameters often change. The vortex dynamics (i.e. the modes and behaviour) influence the stress field and, thus, atherogenic formation. Asymmetry in the flow field determines which wall is more susceptible to atheroma. Similarly, the location of the artery and the degree of stenosis determine the amplitude of flow undulation. Hence, we have conducted the study encompassing a wide range of Reynolds numbers, Womersley numbers and amplitude of the inflow pulsation (described in detail in § 2).
2. Problem description
The geometry consists of an inlet channel of unit non-dimensional height (
$d=1$), which expands to a channel of height, 
$D={{ER}}$, and length, 
$L$, as shown in figure 1. Here 
$L$ is chosen to be sufficiently large to satisfy domain independence as discussed in § 4.1. The figure also represents the notations used for the length of recirculation regions, 
${{X_{r}}}$, which is the axial location along the walls where shear stress is zero. The expansion ratio, 
${{ER}}$, is the only geometrical parameter influencing the nature of the flow.
Schematic of the problem: solid orange line inflow, solid green line outflow, solid black line bounding walls. The shaded area represents the computational domain. Here, 
${{X_{r}}}$ is the length of the recirculation region, where subscripts 1 and 2 represent the bottom and top walls, respectively. The figure denotes an asymmetric flow state, and a Poiseuille profile is prescribed at the inlet.
We consider both steady and pulsatile incompressible flows of a Newtonian fluid through the expansion channel. At the inlet, a sinusoidally varying Poiseuille profile is prescribed of the form
\begin{equation} u(y,t)= \frac{3}{2} (1-4y^{2}) (1 + a \sin (2 {\rm \pi}t /T) ) \quad -1/2 < y < 1/2, \ t \ge 0, \end{equation}
where 
$a$ and 
$T$ are the amplitude and time period of pulsation, respectively.
We define the cross-sectional and temporal average velocity, 
$u_m(t)$ and 
$\bar {u}_m$, respectively as
\begin{equation} u_m(t) = \frac{1}{d} \int_{{-}1/2}^{1/2} u(y,t) \,\mathrm{d}y, \quad \bar{u}_m = \frac{1}{T} \int_{0}^{T} u_m(t) \,\mathrm{d}t. \end{equation}
The inlet profile (2.1) gives a cycle-averaged inflow of unity (
$\bar {u}_m=1$). This leads to the following important hydrodynamic parameters – the reduced velocity 
$U_{r}$, the amplitude 
$a$ and the Reynolds number Re, respectively, given by
\begin{equation} {{U_{r}}}=\bar{u}_m T/d = T, \quad a=\text{amplitude},\quad {\textit{Re}}= \bar{u}_m d/ \nu, \end{equation}
where 
$\nu$ is the kinematic viscosity of the fluid. In the case of steady inflow, the Poiseuille profile 
$u(y)=({3}/{2})(1-4y^{2})$ is directly prescribed. In this case, the Reynolds number, Re, is defined based on the sectional and time-averaged velocity at the inlet as 
${\textit {Re}}= \bar {u}_m d / \nu$.
In this work we examine the stability of pulsating flows through sudden expansion channels for a wide range of 
${{U_{r}}}$, 
$a$ and 
${{ER}}$ (
${{U_{r}}} = 0.5, 1, 1.25, 2.5, 5, 10, 20, 40, 80, 160$; 
${{a}} = 0.25, 0.5, 0.75, 1$; and 
${{ER}} = 2, 3$).
Although we have considered a modulated Poiseuille velocity profile prescribed directly at the expansion inlet, we have also simulated flow through a long inlet channel upstream of the expansion and found that the different modes are similar to that obtained when prescribing the inlet boundary conditions directly at the plane of expansion (these modes are described in § 4.2). However, the critical 
$Re$ values are smaller in long inlet channel cases. For the details, readers are referred to Appendix A. We have also simulated flow using the exact Womersley profile at the plane of expansion as an inlet condition and found similar results in terms of the dynamics of the modes as with time-modulated Poiseuille profile as an inlet at the expansion. For more details, readers are referred to Appendix B.
In the domain of biological flows, a commonly used non-dimensional parameter to denote the flow pulsation (period) is the Womersley number, 
$\alpha$, which is given as 
$\alpha =d/\sqrt {\nu T}$, representing the ratio of the transient inertial force to the viscous force. In this work, we, however, use the reduced velocity 
$U_r$ as a measure of the pulsation time period, and both the quantities are related by 
$\alpha = \sqrt {{\textit {Re}}/ U_r}$. Critical Reynolds numbers obtained in the present study lie between 75 and 400 for 
$ER=2$ for various 
$a$ and 
$U_r$, and correspondingly, Womersley numbers lie between 1.25 and 18.
3. Mathematical modelling
The equations governing the flow are the dimensionless incompressible Navier–Stokes and continuity equation for a Newtonian fluid,
\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}={-}\boldsymbol{\nabla}{p} + \frac{1}{{\textit{Re}}}\Delta \boldsymbol{u}, \end{gather}
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
where 
$\boldsymbol {u}$ consists of the 
$x$ and 
$y$ direction velocity components (
$u$, 
$v$), 
$p$ is the pressure and 
$t$ is the time. The 
$x$ component of the inlet velocity boundary condition is the same as (2.1), with the 
$y$ component velocity being 
$v=0$. The conditions at the outlet are that of a fully developed flow,
\begin{equation} \frac{\partial{\boldsymbol{u}}}{\partial{x}}=0, \end{equation}and, at bounding walls, no-slip conditions are imposed,
\begin{equation} \boldsymbol{u}=0. \end{equation}The equations are solved using a fractional step method (Harlow & Welch Reference Harlow and Welch1965) in a staggered grid arrangement. The convection terms are discretized using a third-order accurate finite difference upwinding scheme, while the diffusion terms use fourth-order central differencing. The time marching is done by Adam–Bashforth's explicit scheme, which is of second-order accuracy. The time step is kept constant at 0.001 for all the simulations, which is found to be lower than 0.1 Courant–Friedrichs–Lewy criteria. The in-house code is accelerated with directive-based programming, OpenACC, to harness the execution capacity on Graphics Processing Units (GPUs). We have used Nvidia V100 16 GB GPU and Intel Xeon Gold 6148 CPU @2.40 GHz in the present work. The linear system of equations resulting from Poisson's equation of pressure is solved by a preconditioned BiCGStab method implemented in the Amgx library (Naumov et al. Reference Naumov2015).
4. Results and discussion
First, we present grid independence and validation study for the steady inlet sudden expansion flow. Validation is done for both symmetric and asymmetric flows.
4.1. Grid independence study and validation
Extensive grid resolution studies have been performed, and the accuracy of the solution is evaluated through Richardson's extrapolation (Slater Reference Slater2021). Richardson's extrapolation is a method to determine the value of a quantity at a continuum limit (grid size being infinitesimal) based on the values at a few discrete grid spacing levels. Order of convergence, 
$P$, deals with the rate at which error diminishes with grid size. In a first-order convergence the solution error scales linearly with the grid size (
$P=1$), while in second-order convergence a reduction in grid size reduces the error quadratically (
$P=2$). Finally, the asymptotic range of convergence, 
$C$, refers to convergence behaviour when grid size is sufficiently small such that changes in grid levels do not affect the solution substantially (
$C=1$). For the definitions of 
$P$, 
$REV$ and 
$C$ please refer to Appendix C.
Table 1 reports the order of convergence, 
$P$, the predicted value of the recirculation length, 
$REV$, using Richardson's extrapolation, and the asymptotic range of convergence, 
$C$. The analysis presented is for three distinct grids with a constant grid ratio, 
$r=2$, for steady inlet at 
${\textit {Re}}=100$ and 
${{ER}}=3$ (see 
$X_{r2}$ versus grid sizes for 
$ER=3$ in table 2).
Grid convergence study. Here 
$P$, 
$REV$ and 
$C$ are the convergence rate, the extrapolated length of the recirculation region and the asymptotic rate of convergence.
Mesh resolution study performed for 
${{ER}}=3$ at 
${\textit {Re}}=100$ and 
${{ER}}=2$ at 
$Re=233.33$.
Table 2 shows the length of the recirculation eddies, 
${{X_{r}}}$, varying with grid sizes for steady Poiseuille inflow at 
${{ER}}=3$ for 
$Re=100$ and 
${{ER}}=2$ for 
$Re=233.33$. In both cases, the Reynolds number was chosen sufficiently above 
${{Re_{c}}}$ (shown later). As shown in the table, a cell size of 
$d/40$ is sufficient enough to resolve the eddies for both 
${{ER}}$ values, and hence, a uniform grid spacing of 
$d/40$ is employed in the present study. Flow features do not change for varying lengths of the domain from 
$L=50$ to 
$L=400$. Hence, the length of the domain is kept sufficiently long, 
$L=100$, for cases with low-to-moderate 
${{U_{r}}}$ (0.5, 1, 1.25, 2.5, 5, 10, 20, 40). For higher 
${{U_{r}}}$ (80, 160, 320), however, the length of the channel is kept at 
$2.5{{U_{r}}}$ to ensure features arising from at least two complete flow cycles reside in the computational domain at any point in time.
For validation, we compare the length of the recirculation region obtained from our simulations using a steady inlet against that of Debuysschère et al. (Reference Debuysschère, Siconolfi, Rimez, Gallaire and Scheid2021). As presented in table 3, the results closely match. We also present the critical Reynolds number obtained for both geometric configurations in table 4, which is in good agreement with the direct numerical simulation (DNS) and LSA results of Debuysschère et al. (Reference Debuysschère, Siconolfi, Rimez, Gallaire and Scheid2021). Furthermore, in figure 2 we show the comparison of the velocity profiles for the case of 
${\textit {Re}}=167$ and 
${\textit{Re}}=233.33$ at 
${{ER}}=2$ with those obtained through Nektar++ (Cantwell et al. Reference Cantwell2015), which is an open-source spectral element method based flow solver. An excellent match is obtained with Nektar++'s prediction. Note the symmetric nature of flow with respect to the 
$y$ axis in figure 2(a) and asymmetric nature in figure 2(b). This is due to the occurrence of a symmetry-breaking bifurcation at 
${{Re_{c}}}=169.5 \pm 0.5$.
Comparison of the length of the recirculation region 
${{X_{r}}}$ at 
${{ER}}=3$ with that of Debuysschère et al. (Reference Debuysschère, Siconolfi, Rimez, Gallaire and Scheid2021).
Critical Reynolds number (
${{Re_{c}}}$) for different expansion ratios (
${{ER}}$). Comparison with DNS and LSA results of Debuysschère et al. (Reference Debuysschère, Siconolfi, Rimez, Gallaire and Scheid2021).
Velocity profiles for (a) 
$Re=167$, (b) 
$Re=233.33$ at 
$ER=2$. At 
$y=0$ from right to left, the velocity profiles are for axial distances 
$x=(1, 2, 4, 6, 10, 25)$. Solid lines represent profiles obtained by our in-house CFD code, while symbols represent those of Nektar++.
In the next section we first introduce the four different kinds of flow modes (denoted as modes 
$A$, 
$B$, 
$C$ and 
$D$) observed in pulsatile sudden expansion channels. Then, we present the mapping of the modes with varying 
${{U_{r}}}$ and 
${{a}}$ and discuss their influence on the critical Reynolds number. The effect of amplitude on the modes is shown next for different values of 
${{U_{r}}}$. Then, we present the mode mapping at a higher 
${{ER}}$. Finally, we explain the modes via the levels of convection and diffusion terms in the vorticity transport equation as well as vorticity production at channel walls and undergo a Floquet analysis to verify the symmetry-breaking 
${{Re_{c}}}$ from our simulations as well as delve into the perturbation flow features responsible for the symmetry breaking.
4.2. General flow features for pulsatile inflow
Here, we report the pulsatile flow features for 
${{ER}}=2$ over a wide range of 
${{U_{r}}}$ (0.5–160) at a constant 
${{a}}$ of 0.25. Then, we investigate various flow modes over a range of inflow amplitudes and at a different expansion ratio.
4.2.1. 
${{U_{r}}}=160, 80, 40$, mode A – synchronized growth of recirculation regions
Figure 3 shows the phase-resolved streamlines at 
${{U_{r}}}=160$, 
$a=0.25$ and 
${\textit{Re}}=170$. At such high reduced velocities (pulsation time period) the eddies behind the step show temporal synchronization with the inflow profile. We see that in the early accelerating part of the cycle, these eddies grow in length (figure 3a,b) and then diminish in the later part of the cycle (figure 3b–d), due to local deceleration at the inlet. We denote this conformal behaviour of the eddies as mode A.
Streamlines for 
${{U_{r}}}=160$ at 
${{a}}=0.25$, and 
${\textit{Re}}=170$ at different phases: (a–d) in steps of 
$T/4$. Here, 
${\textit{Re}}$ is below the critical value and the flow is symmetric. The flow is described by mode A, i.e. conformal growth of recirculation regions with incoming pulsation.
For quantification purposes, we define a scaled recirculation length, 
${{X_{r}^{*}}}$, as 
${{X_{r}^{*}}}={{X_{r}}}/{{\bar {X}_{r}}}$, where 
${{\bar {X}_{r}}}$ denotes the time-average value of 
${{X_{r}}}$. Figure 4(a) shows the temporal variation of 
${{X_{r}^{*}}}$ and the sectional averaged inflow, 
$u_m(t)$, at 
${\textit{Re}}=170$ (
${{U_{r}}}=160$ and 
${{a}}=0.25$). We observe a similar response of the recirculation eddies compared with variation in 
$u_m(t)$ (25 % growth and decay in size about the mean, similar to the amplitude of incoming pulsation). As Re is increased, the flow becomes unstable and the recirculation regions become unequal in length at the top and bottom walls. This asymmetric behaviour (
${{U_{r}}}=160$ and 
${{a}}=0.25$) is shown in figure 4(b) for 
${\textit{Re}}=180$. As plotted against time, the shorter eddy, with a mean length of 4.67 (
${{X}_{{r}2}^{*}}$), becomes phase locked with the inflow profile, while the bigger, lower eddy (with a mean length of 7.1, 
${{X}_{{r}1}^{*}}$) shows an increased phase lag. The phase-locking phenomena of the shorter eddy can be addressed due to its core's close proximity to the inlet since, in this case, the jet first strikes the upper wall before deflecting towards the lower wall. Based on the average flow field, the critical Reynolds number for symmetry-breaking bifurcation, 
${{Re_{c}}}$, is 
$175.5 \pm 0.5$ (
${{U_{r}}}=160$ and 
${{a}}=0.25$), which is slightly higher than that of steady inflow at the same 
${{ER}}$ (
${{Re_{c}}}=169.5 \pm 0.5$).
Scaled recirculation length 
${{X_{r}^{*}}}$ at different phases for 
${{U_{r}}}=160$, 
${{a}}=0.25$ and (a) 
${\textit{Re}}=170$, (b) 
${\textit{Re}}=180$. Results are shown for (c) 
${{X_{r}^{*}}}$ for 
${{U_{r}}}=640, 160, 80, 40$ at 
${{a}}=0.25$ and 
${\textit{Re}}=170$; (d) 
${{X_{r}^{*}}}$ for 
${{U_{r}}}=640, 160, 80, 40$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$. Flow mode A: conformal growth of recirculation regions. Comparison with sectional averaged inflow, 
$u_m(t)$. Plots are dimensionless and shown on the same y-axis scale.
To further explore mode A (at 
${{a}}=0.25$), we observe the response of recirculation regions for other values of 
${{U_{r}}}$ (at 
${\textit{Re}}=170$). When 
${{U_{r}}}$ is large (e.g. 
${{U_{r}}}=640$), the growth of the recirculation region shows negligible phase lag with the inflow profile (see figure 4(c), 
${{X_{r}^{*}}}$ for 
${{U_{r}}}=640$). At 
${{U_{r}}}=160$, the recirculation zone reaches its peak at 
$t/T=5/8$, showing a phase lag of 
$1/8{T}$ with the inflow pulsation. This lag increases to 
$1/4T$ and 
$3/8T$ for 
${{U_{r}}}=80$ and 
${{U_{r}}}=40$, respectively. Note that the recirculation region diminishes in size in the last 
$1/8$th part of the pulsation cycle, whereas the inflow shows an acceleration for all cases, except for the lowest period case (
${{U_{r}}}=640$), in which the response of the recirculation region is nearly synchronized with the inflow. The 
${{Re_{c}}}$ in this regime increases slightly with a reduction in 
${{U_{r}}}$, as 
${{Re_{c}}}=176.5 \pm 0.5$ for 
${{U_{r}}}=80$, and 
${{Re_{c}}}=178.5 \pm 0.5$ for 
${{U_{r}}}=40$ at constant 
${{a}}=0.25$.
The phase lag observed in this case (figure 4c) is similar to the phase lag prevalent in pulsatile flows in a straight channel between the driving pressure gradient and its associated flow rate. As 
$U_r$ decreases, the associated phase lag increases. We define a scaled pressure drop parameter, 
$dp^*(t)=p(10,0,t)-p(0,0,t)/(p(10,0)$
$-\,p(0,0))_{{time\text {-}averaged\ field}}$, where 
$p(10,0,t)$ is considered sufficiently downstream from the end of the recirculation regions. Figure 5 shows the time series plot of 
$dp^*$ for 
$U_r=640, 160, 80, 40$ at 
$a=0.25$ and 
$Re=170$. At 
$U_r=640$, the 
$dp^*$ shows a negligible phase lag, but as 
$U_r$ decreases, the associated phase lag increases. This may result in the phase lag of the growth of the recirculation region as shown in figure 4(c). It can also be observed in figure 5 that as 
$U_r$ decreases, the amplitude of oscillation of 
$dp^*$ increases, suggesting an increased growth and shrinkage of the recirculation region from their respective mean lengths (as observed in figure 4c).
Scaled pressure drop 
$dp^*$ at different phases for 
${{a}}=0.25$ and 
${\textit{Re}}=170$ for various 
${{U_{r}}}=640$, 160, 80 and 40.
4.2.2. ${{U_{r}}}=20$, mode B – necking and diffusion of recirculation region
At 
${{U_{r}}}=20$, 
${{a}}=0.25$, we observe a different flow pattern than discussed in the previous subsection. The streamlines are depicted in figure 6 for 
${\textit{Re}}=180$ (which is below the 
${{Re_{c}}}$ for the chosen pulsation parameters). There is a large growth in the size of the recirculation region during a major part of the inflow pulsation cycle (
$t/T=0 \text { to } 3/4$; see figure 6a–d), and then it diminishes rapidly over a very short duration of time (
$t/T=3/4\text { to }1$; see figure 6d–a).
Streamlines for 
${{U_{r}}}=20$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$ at different phases (a–g). Flow mode B: necking and diffusion of the recirculation region.
This abrupt reduction in the size of the recirculation region occurs via a necking phenomenon (here it occurs at 
$t/T=0.86$; see figure 6), which splits the recirculation region into a proximal (at the corner of the expansion) and distal portion. The distal portion then diffuses rapidly into the flow. The scaled recirculation length is shown in figure 7, which signifies the abrupt decline in the size of the recirculation region after reaching maximum growth. This necking and diffusion of the recirculation zone is termed mode B here. The critical Reynolds number in this regime is higher than that of mode A, i.e. for 
${{U_{r}}}=20$ and 
${{a}}=0.25$, 
${{Re_{c}}}=187.5 \pm 0.5$. We also plot the scaled pressure drop, 
$dp^*$, in this case (see figure 7b). The scaled pressure drop changes its sign from being adverse (
$dp^*=+{\rm ve}$) to being favourable (
$dp^*=-{\rm ve}$) at 
$t/T=0.8$, during the lower inflow velocity phase, which corresponds to the necking behaviour of the recirculation zones at 
$t/T=0.86$.
(a) Scaled recirculation length 
${{X_{r}^{*}}}$ at different phases for 
${{U_{r}}}=20$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$. Comparison with sectional averaged inflow, 
$\bar {u}_m(t)$. Parameters, 
$\bar {u}_m(t)$ and 
${{X_{r}^{*}}}$, are dimensionless and shown on the same axis scale. (b) Scaled pressure drop 
$dp^*$ at different phases for 
${{a}}=0.25$ and 
${\textit{Re}}=180$ for 
${{U_{r}}}=20$.
Above 
${{Re_{c}}}$, asymmetricity in this regime is observed throughout the cycle. The necking takes place at two locations in the eddy associated with the upper wall while a single necking occurs in that of the lower wall beyond 
${{Re_{c}}}$. This necking phenomenon can be attributed to an increased convection and smaller diffusion rates of vorticity. This will be discussed in a later subsection.
4.2.3. ${{U_{r}}}=10, 5$, mode C – splitting of recirculation region
In this high-frequency pulsatile flow regime there is a similar necking phenomenon as discussed in mode B; however, the distal eddy advects downstream instead of diffusing in the flow in this mode. Flow behaviour for this case is shown as a streamline plot in figure 8 for 
${{U_{r}}}=10$. The recirculation eddies generated at the corner of expansion advects to a distal position after necking and subsequent breakdown (see 
$x=5.5$ in figure 8(a); the necking occurs in the last 
$1/8$th part of the cycle as in the previous subsection). The advected eddy diffuses later through the subsequent cycles. We do not provide the length of the recirculation zone here because the recirculation eddies split and convect. Rather, we present the scaled location of the centre of the eddy 
${{X_{c}^{*}}}$, (
${{X_{c}^{*}}}=X_c/\bar {X}_{c}$, where 
$\bar {X}_c$ denotes the time-average value of 
$X_c$) in figure 9(a). We observe that the core position of the recirculation zone moves at a constant speed. The scaled pressure drop, 
$dp^*$, shown in figure 9(b), shows a higher amplitude compared with the two low-frequency modes (A and B) discussed earlier. The critical Reynolds number for 
${{U_{r}}}=10$ and 
${{U_{r}}}=5$ is 
${{Re_{c}}}=201.5 \pm 0.5$ and 
${{Re_{c}}}=180.5 \pm 0.5$, respectively. This splitting of the recirculation region is termed mode C for subsequent discussion. A similar observation of the recirculation region splitting and convection is found in pulsatile flow over a constricted tube (figure 2 in Barrere et al. Reference Barrere, Brum, Anzibar, Rinderknecht, Sarasua and Cabeza2023). Higher rates of vortex convection in this regime are attributed to both mean and fluctuating components of convection of vorticity (discussed in a later subsection).
Streamlines for 
${{U_{r}}}=10$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$ at different phases: (a–d) phases in steps of 
$T/4$. Flow mode C: splitting and advection of the recirculation region.
(a) Scaled recirculation region centre position 
$X_c^*$ at different phases for 
${{U_{r}}}=10$, 
${\textit{Re}}=180$ and 
${{a}}=0.25$. Comparison with sectional averaged inflow, 
$\bar {u}_m(t)$. Parameters, 
$\bar {u}_m(t)$ and 
$X_c^*$, are dimensionless and shown on the same axis scale. (b) Scaled pressure drop 
$dp^*$ at different phases for 
${{a}}=0.25$ and 
${\textit{Re}}=180$ for 
${{U_{r}}}=10$.
4.2.4. ${{U_{r}}}=2.5, 1.25, 1, 0.5$, mode D – inverse growth of recirculation region
When 
${{U_{r}}}$ is small, we observe a high phase lag in the behaviour of eddies with respect to the inflow pulsation. Figure 10 shows the streamlines in this regime for 
${{U_{r}}}=1.25$ and 
${\textit{Re}}=160$ at 
${{a}}=0.25$. During acceleration (figures 10a,b and 10d,a), the size of the recirculation region shrinks, while during deceleration (figure 10b–d) it grows, showing an anti-phase behaviour with inflow pulsation. The critical Reynolds number for this case is 
$161.5 \pm 0.5$, which is lower than the steady inflow case at this 
${{ER}}$. The temporal variation of 
${{X_{r}^{*}}}$ for different 
${{U_{r}}}$ is shown in figure 11(a), and it shows a complete off-phase nature as compared with the variation of sectional averaged inflow. The scaled pressure drop, 
$dp^*$, for the case of 
${{U_{r}}}=1.25$, is shown in figure 11(b). We can observe a very high amplitude of the pressure drop compared with the previous cases. When the pressure drop is favourable (
$dp^*=-{\rm ve}$, 
$0 \le t/T \le 0.25$), the recirculation region shrinks in size (see figure 11(a), 
$0 \le t/T \le 0.25$). Then, when the pressure drop is adverse (
$dp^*=+{\rm ve}$, 
$0.25 \le t/T \le 0.75$), the recirculation zone grows in size. Finally, the pressure drop becomes favourable (
$dp^*=+{\rm ve}$, 
$0.75 \le t/T \le 1$), resulting in a decline in the size of the recirculation region. Therefore, the recirculation zone shrinks and grows following the sign of the pressure drop and, hence, shows phase lag with both inflow pulsation and pressure drop variation cycles. This subsection describes the fourth mode, i.e. mode D – showing inverse growth of the recirculation region. Compared with other modes, this mode shows a higher diffusion and comparable convection rates of vorticity. We have also observed higher rates of convection of fluctuating vorticity (described in § 4.6).
Streamlines for 
${{U_{r}}}=1.25$ at 
${{a}}=0.25$ and 
${\textit{Re}}=160$ at different phases: (a–d) phases in steps of 
$T/4$. Flow mode D: inverse growth of the recirculation region.
(a) Scaled recirculation length 
${{X_{r}^{*}}}$ at different phases for 
${{U_{r}}}=2.5, 1.25, 1, 0.5$ at 
${{a}}=0.25$ and 
${\textit{Re}}=160$. Comparison with sectional averaged inflow, 
$\bar {u}_m(t)$. Parameters, 
$\bar {u}_m(t)$ and 
${{X_{r}^{*}}}$, are dimensionless and shown on the same axis scale. (b) Scaled pressure drop 
$dp^*$ at different phases for 
${{a}}=0.25$ and 
${\textit{Re}}=160$ for 
${{U_{r}}}=1.25$.
The above results describe in detail the nature of various flow modes at 
${{a}}=0.25$. These modes encompass the flow description at other 
${{a}}$ and 
${{ER}}$ as well. It is to be mentioned that the modes for different 
$a$ and 
$U_r$ regimes are similar for a long inlet channel flow as well as for a Womersley inlet velocity profile as discussed in Appendices A and B.
4.3. Critical Reynolds number and mode map
Figure 12(a) shows the 
${{Re_{c}}}$ associated with differing 
${{U_{r}}}$ and 
${{a}}$ along with the modes. At low 
${{U_{r}}}$ 
$(0.5, 1, 1.25, 2.5)$, mode D (inverse growth of recirculation region) is prevalent, with increasing 
${{a}}$ leading to a reduction in stability. At moderate 
${{U_{r}}}$ 
$(5, 10, 20)$, mode C (splitting and advection of the recirculation region) is seen, which is associated with higher stability. At high 
${{U_{r}}}$ 
$(40, 80, 160)$, mode B (necking and diffusion of the recirculation region) is prevalent with some of the highest stability points, and mode A (conformal growth of the recirculation region) is also observed to stabilize the flow at further higher 
${{U_{r}}}$. At low or moderate reduced velocities (
${{U_{r}}} \leq 5$), stability decreases with 
${{a}}$. Above 
${{U_{r}}}=10$, stability increases with 
${{a}}$ (0.25, 0.50, 0.75) except at 
${{a}}=1.00$, which is reported to have lower stability than 
${{a}}=0.75$. This may be attributed to the presence of a zero inflow phase. The mode map is depicted separately in figure 12(b), which shows that a higher time period (or 
${{U_{r}}}$) yields mode A at low or moderate amplitudes. However, mode B is observed at high values of 
${{a}}$. As 
${{U_{r}}}$ reduces, mode A is only observed in lower 
${{a}}$, and moderately higher 
${{a}}$ yields mode B. As 
${{U_{r}}}$ further decreases, mode C is also seen at high 
${{a}}$. At much lower 
${{U_{r}}}$, we see only mode C across all values of 
${{a}}$. As 
${{U_{r}}}$ is lower than 
$5$, only mode D is seen. We have not investigated beyond 
${{a}}=1.00$ to avoid backflow.
(a) Critical Reynolds number 
${{Re_{c}}}$, annotated with flow modes, (b) mode map, for 
${{U_{r}}}$ (
$0.5\unicode{x2013}160$) and 
${{a}}$ (
$0.25$–
$1.00$) for 
$ER=2$.
4.4. Effect of amplitude, $a$, on the flow features
4.4.1. Low ${{U_{r}}}$ $(0.5, 1, 1.25, 2.5)$
At low 
${{U_{r}}}$, as 
${{a}}$ increases, the disappearance of the recirculation region speeds up in mode D. Figure 13(a) shows an increase in the opposite response of 
${{X_{r}^{*}}}$ with sectional averaged inflow, 
$u_m(t)$, at increased 
${{a}}$. The dashed lines represent the growth of the recirculation regions until it spans over the entire domain length at higher values of 
${{a}}$, while filled circles represent the reoccurrence of the recirculation at later times. This particular flow field is visualized in figure 14 for 
${{a}}=0.75$ at 
${{U_{r}}}=1$ and 
${\textit{Re}}=108$. The critical Reynolds number decreases from 162.5 to 86.5 when 
${{a}}$ is increased from 0.25 to 1.00. Here, 
${{Re_{c}}}$ for the steady case is 169.5.
(a) Plot of 
${{X_{r}^{*}}}$ at different phases for 
${{U_{r}}}=1$, for various 
${{a}}$ and 
$Re$ just below respective 
${{Re_{c}}}$. The magnitude of the inverse response increases with 
${{a}}$. Sharp changes of curvature of 
${{X_{r}^{*}}}$ at the deceleration part for 
${{a}}$ (0.50, 0.75 and 1.00) are indicative of back feeding of the recirculation zone from the outlet section. (b) Plot of 
${{X_{r}^{*}}}$ at different phases for 
${{U_{r}}}=160$, for 
${{a}}=0.25, 0.50, 0.75, 1.00$ and Re is just below respective 
${{Re_{c}}}$ (
${\textit{Re}}=170$ for 
${{a}}=0.25$; 
${\textit{Re}}=185$ for 
${{a}}=0.50$; 
${\textit{Re}}=228$ for 
${{a}}=0.75$; 
${\textit{Re}}=269$ for 
${{a}}=1.00$). The coloured lines show the flow mode changes from A to B with increasing a.
Streamlines for 
${{U_{r}}}=1$ at 
${{a}}=0.75$ and 
${\textit{Re}}=108$ at different phases: (a–c) phases in steps of 
$T/8$, (a) representing recirculation regions spanning over the entire domain length, (b) representing the shrinkage of the recirculation region in the wall normal direction and (c) representing reoccurrence of the same at a later time.
4.4.2. High ${{U_{r}}}$
When 
${{U_{r}}}$ is high (
${{U_{r}}}=160$), higher values of 
${{a}}$ augment the stability of the flow. This happens due to a change of modes from A to B as 
${{a}}$ is increased (see figure 13b). For 
${{a}} = 0.25$, 0.50, and 0.75, mode A with conformal growth and shrinkage of the recirculation region is prevalent. On the other hand, at 
${{a}}=1.00$, the mode is necking and diffusion of the recirculation region (mode B), with a sharp reduction in the size of the recirculation region. This is associated with improved stability at 
${{a}}= 1.00$ at 
${{U_{r}}}=160$. Here 
${{Re_{c}}}$ increases from 175.5 to 245.5 as 
${{a}}$ increases at this 
${{U_{r}}}$.
4.4.3. Moderate ${{U_{r}}}$
At moderate values of 
${{U_{r}}}$ (5–20), mode C is prevalent. Interestingly, at 
${{U_{r}}}=10$, there is a change of trend, whereby, for lower values of 
${{U_{r}}}$, the 
${{Re_{c}}}$ decreases with 
${{a}}$, but at higher 
${{U_{r}}}$, the same feature is no longer maintained.
4.5. Critical Reynolds number and mode map at a higher ${{ER}}$
Now, we report the mode map and critical Reynolds number for varying 
${{U_{r}}}$ and 
${{a}}$ at 
${{ER}}=3$. Figure 15(a,b) shows the critical Reynolds numbers and associated mode map, respectively. A similar trend is seen in both the 
${{ER}}$ 
$(2, 3)$. The 
${{Re_{c}}}$ are lower than those of 
${{ER}}=2$, and the associated mode map shows a slight change in some positions, especially in moderate 
${{U_{r}}}$. Mode A is delayed to high 
${{U_{r}}}$ here. The reduction in the values of the critical Reynolds number at a higher expansion ratio is reported for steady inflow cases (Battaglia et al. Reference Battaglia, Tavener, Kulkarni and Merkle1997; Drikakis Reference Drikakis1997). The same effect may persist in the pulsatile cases, lowering the 
${{Re_{c}}}$.
(a) Critical Reynolds number 
${{Re_{c}}}$, annotated with flow modes, and (b) mode map, for 
${{U_{r}}}$ (
$0.5\unicode{x2013}160$) and 
${{a}}$ (
$0.25$–
$1.00$) for 
${{ER}}=3$.
4.6. Transport of vorticity in different modes
In earlier subsections we have described in detail the various flow modes in sudden expansion channels with inlet pulsatility depending on the amplitude and frequency of the inflow waveform. In this section we explain the vorticity dynamics associated with the respective modes. Consider the vorticity transport equation
\begin{equation} \frac{\partial{\omega}}{\partial{t}} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\omega = \frac{1}{{\textit{Re}}}\Delta{\omega}, \end{equation}
where 
$\omega$ denotes the vorticity in the 
$z$ direction, the second term on the left-hand side of (4.1) denotes convection of vorticity, and the right-hand side term denotes diffusion of vorticity. We have also considered the production of vorticity at the channel walls for all the modes. For that purpose, vorticity flux at the channel wall is estimated from the pressure gradient along the channel wall as
\begin{equation} \frac{\partial p}{ \partial x}={-}\mu \frac{\partial \omega}{\partial y}. \end{equation}
It may be inferred from (4.2) that the positive value of the pressure gradient along the lower channel wall (
${\partial p}/{\partial x} > 0$) indicates the growth of a negative shear layer over that wall, and its negative value indicates weakening of this shear layer. For the top wall, the wall normal being in the negative 
$y$ direction, the effects are similar for the positive shear layer formed over that wall. We have discussed the variation of the wall pressure gradients during different phases of pulsation and its effects on the respective modes in the following subsections.
4.6.1. Mode A
Figures 16(a), 16(b) and 16(c) show the contours of instantaneous vorticity, instantaneous convection and diffusion of vorticity for mode A, respectively, at 
$t/T=0$, 
$t/T=0.125$ and 
$t/T=0.71875$ for 
${{a}}=0.75$, 
${{U_{r}}}=160$ and 
${\textit{Re}}=228$. We observe an enhanced absolute convection of vorticity relative to the diffusion term (see figure 16b) during the high inlet velocity accelerating phase. Therefore, (4.1) provides an increase in the rate of change of vorticity (positive in the upper shear layer and negative in the lower shear layer), which elongates and strengthens the shear layers and, consequently, results in an increase in the size of the recirculation region. Later, as the inlet velocity reduces, the convection term shows a diminished magnitude, but diffusion increases due to higher gradients of vorticity, which was established by the shear layers. We also observe enhanced diffusion of vorticity into walls in the later phases (see figure 16c). Finally, the diffusion of vorticity becomes greater than convection, which results in a net reduction rate of change of vorticity, and the shear layer weakens along with shrinkage in the size of the recirculation region (see figure 16c). In this mode the vorticity convection term shows a smooth spatial distribution, indicating no break-up of the recirculation region. We observe synchronized growth and decay of the recirculation region.
Contours for vorticity, convection and diffusion of vorticity for mode A for 
${{U_{r}}}=160$, 
${{a}}=0.75$ and 
$Re=228$ at (a) 
$t/T=0$, (b) 
$t/T=0.125$ and (c) 
$t/T=0.71875$ phase.
Figure 20(a) shows the pressure gradient along the channel walls from which vorticity flux through the wall may also be estimated (refer to (4.2)). As the recirculation zone stays over the wall during the entire pulsation cycle without showing necking or disappearance, a positive value of wall pressure gradient is observed in this region over all phases. However, we can observe that for mode A, 
$\partial p/ \partial x$ increases from 
$t/T=0$ to 
$t/T=0.25$ and then it decreases till 
$t/T=0.75$ and then it again increases to 
$t/T=1$ (figure 20a), giving rise to a similar nature of rate of change of vorticity influx from the wall and, consequently, the conformal nature of growth and decay of the recirculation regions (growth during the first quarter, then decay and again showing growth at the last quarter) is observed in mode A.
4.6.2. Mode B
Mode B (see figure 17a–c) is typically interesting as it encompasses the phenomena of eddy splitting. Here, at a later part of the pulsatile cycle (
$t/T=0.8325$), there is a discontinuity in the contours of convection of vorticity with its peaks at two locations (one at the vortex tip, another at the lip of expansion, see figure 17b,c). Due to this patchy convection term, the shear layer does not grow smoothly. Furthermore, higher vorticity diffusion is observed between the peaks of the convection term that leads to the necking and eventual splitting of the recirculation region in this mode. Figure 17(c) shows the necking in the recirculation bubble.
Contours for vorticity, convection and diffusion of vorticity for mode B for 
${{U_{r}}}=20$, 
${{a}}=0.25$ and 
$Re=180$ at (a) 
$t/T=0$, (b) 
$t/T=0.655$, and (c) 
$t/T=0.8325$ phases.
For mode B, the pressure gradient along the wall remains positive for a longer span in the axial direction at 
$t/T=0.50$, indicating an increased span of the recirculation region till 
$t/T=0.50$ (figure 20b). At 
$t/T=0.75$, however, the pressure gradient falls to zero value after a short span and again rises to the second peak, which is higher than the first. This behaviour indicates a local reduction in vorticity influx within the recirculation region, representing the necking phenomena in this mode.
4.6.3. Mode C
In mode C (see figure 18a) the eddy-splitting phenomenon results in the convection of distinct recirculation regions. We can see reduced diffusion in the vicinity of the core of the recirculation region, yielding this distinct vortex feature. To further explore the vorticity dynamics during the convection of the vortices in this high-frequency pulsatile flow, we look into the time-averaged vorticity equation following the decomposing of the flow field in mean and fluctuating components. Considering 
$\boldsymbol {u}=\bar {\boldsymbol {u}} + \hat {\boldsymbol {u}}$ and 
$\omega =\bar {\omega } + \hat {\omega }$, where the instantaneous fields (
$\boldsymbol {u}(x,y,t)$, 
$\omega (x,y,t)$) are decomposed into a mean (
$\bar {\boldsymbol {u}}(x,y)$, 
$\bar {\omega }(x,y)$) and an oscillating field (
$\hat {\boldsymbol {u}}(x,y,t)$, 
$\hat {\omega }(x,y,t)$) for pulsatile flow, and averaging (4.1) over a time period, the averaged vorticity transport equation is
\begin{equation} \frac{\partial{\bar{\omega}}}{\partial{t}} + \underbrace{(\bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla})\bar{\omega}}_{\text{term I}} + \underbrace{\overline{(\hat{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla})\hat{\omega}}}_{\text{term II}} = \underbrace{\frac{1}{{\textit{Re}}}\Delta{\bar{\omega}}}_{\text{term III}}. \end{equation}In the above equation, term I represents the convection of averaged vorticity by the mean flow field, term II represents averaged convection of fluctuating vorticity by the fluctuating flow field and term III represents the diffusion of averaged vorticity. In the time-averaged vorticity transport plot (figure 18b), we note that transport term III is stronger than the other two terms and shows distinct peaks at the tips of the time-averaged recirculation regions. Therefore, the unsteadiness in both vorticity and velocity is responsible for the departure of the recirculation region through this location.
(a) Contours for vorticity, convection and diffusion of vorticity for mode C for 
${{U_{r}}}=20$, 
${{a}}=0.75$ and 
$Re=342$ at 
$t/T=0$; (b) contours of averaged vorticity transport rates for the 
$t/T=0$ phase.
For mode C, the peaks and valleys in the pressure gradient are responsible for necking and advection of recirculation zones (figure 20c).
4.6.4. Mode D
Figures 19(a), 19(b) and 19(c) show the contours of instantaneous vorticity, instantaneous convection and diffusion of vorticity for mode D, respectively, at 
$t/T=0$, 
$t/T=0.25$ and 
$t/T=0.50$, for 
${{U_{r}}}=1$, 
${{a}}=0.75$ and 
${\textit{Re}}=108$. In this mode we have an inverse relation of the growth of the recirculation region with the inflow pulsation. As the inflow flux increases (
$t/T=0$ to 
$t/T=0.25$), the diffusion of vorticity towards the wall is higher (see figure 19a,b). Owing to this diffusion behaviour, the shear layer grows closer to the wall, and consequently, we observe attached boundary layers over the wall with the disappearance of the recirculation bubble during the higher velocity phases of the pulsation cycle. Convection of vorticity shows a patchy anti-symmetric behaviour, and its magnitude is not substantial away from the inlet. Therefore, the shear layer does not spread further downstream, unlike the previous modes. Small local fluctuations in the shear layer are responsible for the patches in the convection term. At the low-velocity phases of the cycle (
$t/T=0.50$ and later), the resultant vorticity from the wall diffuses back to the core shear layer (see figure 19c). This results in the movement of the shear layer away from the wall and the eventual reappearance of the recirculation regions at the deceleration phase of the cycle. Overall, this mode shows an inverse nature of decay and growth of the recirculation region as compared with velocity acceleration and deceleration.
Contours for vorticity, convection and diffusion of vorticity for mode D for 
${{U_{r}}}=1$, 
${{a}}=0.75$ and 
$Re=108$ at (a) 
$t/T=0$, (b) 
$t/T=0.25$ and (c) 
$t/T=0.50$ phases.
Pressure gradient at bottom wall at 
$t/T=0$, 
$t/T=0.25$, 
$t/T=0.50$ and 
$t/T=0.75$ for (a) mode A, (b) mode B, (c) mode C and (d) mode D, respectively.
In mode D, at 
$t/T=0$ (figure 20d), the pressure gradient is negative, signifying a positive vorticity influx to the negative shear layer formed over the bottom wall, and hence, the recirculation region shrinks in span and size from 
$t/T=0$ to 
$t/T=0.25$. The positive shear layer of the top surface shrinks similarly due to the influx of negative vorticity there. From 
$t/T=0.25$ to 
$t/T=0.5$, the pressure gradient is positive and results in an influx of vorticity from the walls, which strengthens the shear layers and yields an increased size of the recirculation region. This behaviour also explains the inverse growth of the recirculation region with respect to the inflow pulsation.
4.7. Floquet analysis
In this subsection we explore the stability of time-periodic symmetric flow in sudden expansion via Floquet analysis. A similar analysis has been performed by Sherwin & Blackburn (Reference Sherwin and Blackburn2005) and Gopalakrishnan, Pier & Biesheuvel (Reference Gopalakrishnan, Pier and Biesheuvel2014) for pulsatile flow over a model stenotic artery and model abdominal aortic aneurysms, respectively.
4.7.1. Perturbation equations
In Floquet analysis one examines the stability of spatially developed and time-periodic base flows over symmetric and antisymmetric perturbations, and one expects the antisymmetric perturbations to become unstable and lead to potential symmetry breaking in our configuration. Let the composite velocity field be decomposed as
\begin{equation} \boldsymbol{u}=\boldsymbol{U}+\boldsymbol{u'}, \end{equation}
where 
$\boldsymbol {U}(x,y,t)$ is the 
$T$-periodic base flow whose stability is to be determined, and 
$\boldsymbol {u'}(x,y,t)$ is an infinitesimal two-dimensional perturbation imposed on the base flow.
The 
$T$-periodic base flow is obtained by solving the respective variables 
$(\boldsymbol {U},P)$ in the Navier–Stokes equation,
\begin{equation} \frac{\partial{\boldsymbol{U}}}{\partial{t}} + (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{U}={-}\boldsymbol{\nabla}{P} + \frac{1}{{\textit{Re}}}\Delta \boldsymbol{U}, \quad \text{with} \ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U=0} \quad \text{in} \ \varOmega/2\ \text{(half of domain)}, \end{equation}
the boundary conditions for 
$\boldsymbol {U}$ are obtained as consistent with that of 
$\boldsymbol {u}$ ((2.1), (3.3), (3.4)).
Substituting (4.5) in (3.1) and linearizing with respect to perturbation terms results in the following equation for the perturbation quantities:
\begin{equation} \frac{\partial{\boldsymbol{u'}}}{\partial{t}} ={-} (\boldsymbol{u'} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{U} - (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u'} -\boldsymbol{\nabla}{p'} + \frac{1}{{\textit{Re}}}\Delta \boldsymbol{u'}, \quad \text{with}\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u'}=0 \quad\text{in} \ \varOmega. \end{equation}
At all boundaries, 
$\boldsymbol {u'}=0$ is specified.
Defining an operator 
$\boldsymbol {L}(\boldsymbol {u'})$ consisting of the right-hand side of the linearized equation (4.6),
\begin{equation} \frac{\partial{\boldsymbol{u'}}}{\partial{t}} = \boldsymbol{L}(\boldsymbol{u'}), \end{equation}
is an evolution equation based on 
$\boldsymbol {u'}$ and is of Floquet type. The operator 
$\boldsymbol {L}(\boldsymbol {u'})$ is also T-periodic, and 
$\boldsymbol {u'}(x,y,t) \equiv \boldsymbol {\tilde {u}}(x,y,t) \exp (\sigma t)$, where 
$\boldsymbol {\tilde {u}}(x,y,t)$ are also T-periodic functions. The Floquet modes are the eigenfunctions of the operator 
$\boldsymbol {L}$, and the complex numbers 
$\sigma$ are the Floquet exponents. The velocity and pressure perturbations from cycle to cycle are related by
\begin{equation} \boldsymbol{u'}(x,y,t+T)=\exp (\sigma T) \boldsymbol{u'}(x,y,t). \end{equation}
Thus, the perturbations may grow or decay exponentially as per the sign of 
$\sigma$, which is also known as (complex) growth rate. Floquet multipliers are related to Floquet exponents by 
$\mu = \exp (\sigma T)$. For 
$|\mu | > 1$, the flow is unstable, and for 
$|\mu | < 1$, the flow is stable. A value of 
$|\mu | = 1$ represents neutral stability. In our work, we came across a critical (real) Floquet multiplier, 
$\mu =+1$, resulting in a synchronous instability with the same period 
$T$. We have neither seen the occurrence of period-doubling bifurcation 
$\mu =-1$ nor Neimark–Sacker bifurcation 
$\mu =\exp (\pm \iota \theta )$.
The symmetric base flow is obtained by simulating half of the channel with symmetric boundary conditions (
${\partial u}/{\partial y}=0$ and 
$v=0$ at 
$y=0$) imposed at the centreline and then reflecting the half-domain results in a complete channel.
Validation of Floquet multiplier 
$\mu$ at 
$Re=190$, 
$\beta =1.6$ with Barkley & Henderson (Reference Barkley and Henderson1996). Here, 
$N_p$ is the polynomial order.
Convergence study of Floquet multiplier 
$\mu$ as a function of polynomial order 
$N_p$ for 
${{U_{r}}}=1$, 
$a=0.25$ and 
$Re=180$. For each 
$N_p$, 
$N_{tot}$ is the number of independent mesh points for the spectral element mesh.
4.7.2. Numerical procedure
The stability analysis is performed by a Krylov-subspace iteration of randomized perturbations through (4.7). An Arnoldi method is used to extract dominant eigenpairs. For all cases, 256 time slices are provided from the T-periodic base flow, and intermediate stages are approximated through Fourier-series reconstruction. The stability analysis is performed in open-source Nektar++ (Cantwell et al. Reference Cantwell2015) software. For more details on the numerical procedure, readers are referred to the earlier works of Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008) and Tuckerman & Barkley (Reference Tuckerman, Barkley, Doedel and Tuckerman2000).
4.7.3. Validation
We have considered the benchmark case of the onset of three-dimensional instability in periodic vortex shedding over a circular cylinder for the validation of the Floquet eigenvalue solver. Table 5 compares the Floquet multiplier value for flow over a cylinder near the transition from two-dimensional periodic flow (von Kármán street) to the flow being unstable to three-dimensional perturbations (
$Re_c=188.5 \pm 1$) with that of Barkley & Henderson (Reference Barkley and Henderson1996) at 
$Re=190$. An excellent match is obtained for polynomial order 8 and above.
Validation of 
${{Re_{c}}}$ from Floquet study and simulation.
4.7.4. Convergence study for the present problem
A series of convergence tests were performed in order to determine the appropriate polynomial order, 
$N_p$, for the base flow and Floquet stability analysis. Table 6 shows the convergence test for the leading Floquet multiplier with polynomial order for 
${{U_{r}}}=1$, 
$a=0.25$ and 
$Re=180$. As shown from 
$N_p=4$ to 
$N_p=8$, there is a deviation of 0.003 %, and hence, we perform the remainder of the stability study using 
$N_p=4$, keeping in view the computational costs.
Critical Reynolds number (
${{Re_{c}}}$) for 
${{ER}}=2$ for a long inlet channel upstream of the expansion region. Comparison with Drikakis (Reference Drikakis1997).
4.7.5. Perturbation flow features
In this subsection we further explore the asymmetricity in flow behaviour beyond the critical Reynolds number in each of the flow modes (A, B, C and D). We also show the velocity perturbations responsible for the asymmetric flow nature via a Floquet analysis along with the respective Floquet multipliers. The presented results are above 
${{Re_{c}}}$, representing asymmetric flow in each of the cases. Figure 21 shows the streamline flow on top of 
$z$ vorticity contours for 
${{U_{r}}}=40, {{a}}=1.25\ \text {and } {\textit{Re}}=180$ at 
$0$th phase, describing the asymmetric nature of flow for mode A. Throughout the flow pulse, the flow remains asymmetric, with the upper eddy remaining longer than the lower one. The figure also accompanies the velocity perturbations 
$\tilde {u}$ and 
$\tilde {v}$ of the most unstable eigenpairs. The eigenfunction 
$\tilde {u}$ is antisymmetric with respect to the centreline, while 
$\tilde {v}$ is symmetric (figure 21). Consequently, we can see the same sign of perturbation vorticity near both top and bottom walls, which results in elongation of the upper shear layer and shrinking of the bottom shear layer, leading to the overall asymmetry. The Floquet multiplier obtained is 1.32 that is greater than 1, representing an unstable symmetric base flow at this Re and gives rise to an asymmetric state. The figure represents unstable mode A.
(a) Streamlines on top of vorticity contours for 
${{U_{r}}}=40$, 
${{a}}=0.25$ and 
$Re=180$ at 
$t/T=0$ for unstable mode A. (b) The 
$x$ and (c) 
$y$ components of the most unstable velocity perturbations and (d) vorticity perturbation. Floquet multiplier, 
$\mu =1.32$.
(a) Streamlines on top of vorticity contours for 
${{U_{r}}}=20$, 
${{a}}=0.25$ and 
$Re=200$ at 
$t/T=0$ for unstable mode B. (b) The 
$x$ and (c) 
$y$ components of the most unstable velocity perturbations and (d) vorticity perturbation. Floquet multiplier, 
$\mu =1.12$.
(a) Streamlines on top of vorticity contours for 
${{U_{r}}}=10$, 
${{a}}=0.25$ and 
$Re=220$ at 
$t/T=0$ for unstable mode C. (b) The 
$x$ and (c) 
$y$ components of the most unstable velocity perturbations and (d) vorticity perturbation. Floquet multiplier, 
$\mu =1.10$.
(a) Streamlines on top of vorticity contours for 
${{U_{r}}}=1$, 
${{a}}=0.25$ and 
$Re=180$ at 
$t/T=0$ for unstable mode D. (b) The 
$x$ and (c) 
$y$ components of the most unstable velocity perturbations and (d) vorticity perturbation. Floquet multiplier, 
$\mu =1.01$.
Streamlines for 
${{U_{r}}}=160$ at 
${{a}}=0.25$ and 
${\textit{Re}}=170$ at different phases: (a–d) in steps of 
$T/4$ with a long channel and modulated Poiseuille imposed at the inlet. The flow is asymmetric and described by mode A, i.e. conformal growth of the recirculation region.
Streamlines for 
${{U_{r}}}=20$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$ at different phases: (a–d) in steps of 
$T/4$ with a long channel and modulated Poiseuille imposed at the inlet. (e,f) Selected time points showing necking phenomena. The flow is asymmetric and described by mode B, i.e. necking and diffusion of the recirculation region.
Streamlines for 
${{U_{r}}}=10$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$ at different phases: (a–d) in steps of 
$T/4$ with a long channel and modulated Poiseuille imposed at the inlet. The flow is asymmetric and described by mode C, i.e. splitting and advection of the recirculation region.
Streamlines for 
${{U_{r}}}=1.25$ at 
${{a}}=0.25$ and 
${\textit{Re}}=160$ at different phases: (a–d) in steps of 
$T/4$ with a long channel and modulated Poiseuille imposed at the inlet. The flow is asymmetric and described by mode D, i.e. inverse growth of the recirculation region.
Comparison of velocity profiles at inlet (plane of expansion) between exact Womersley and modulated Poiseuille for 
$U_r=160, a=0.25, Re=170$. Results are shown for (a) 
$t/T=0$, (b) 
$t/T=0.25$, (c) 
$t/T=0.50$, (d) 
$t/T=0.75$, (e) time averaged.
Streamlines for 
${{U_{r}}}=160$ at 
${{a}}=0.25$ and 
${\textit{Re}}=170$ at different phases: (a–d) in steps of 
$T/4$ with an exact Womersley profile imposed at the inlet. The flow is described by mode A.
Comparison of velocity profiles at inlet (plane of expansion) between exact Womersley and modulated Poiseuille for 
$U_r=20, a=0.25, Re=180$. Results are shown for (a) 
$t/T=0$, (b) 
$t/T=0.25$, (c) 
$t/T=0.50$, (d) 
$t/T=0.75$, (e) time averaged.
Streamlines for 
${{U_{r}}}=20$ at 
${{a}}=0.25$, and 
${\textit{Re}}=180$ at different phases: (a–d) in steps of 
$T/4$ with an exact Womersley profile imposed at the inlet. The flow is described by mode B.
Comparison of velocity profiles at inlet (plane of expansion) between exact Womersley and modulated Poiseuille for 
$U_r=10, a=0.25, Re=180$. Results are shown for (a) 
$t/T=0$, (b) 
$t/T=0.25$, (c) 
$t/T=0.50$, (d) 
$t/T=0.75$, (e) time averaged.
Streamlines for 
${{U_{r}}}=10$ at 
${{a}}=0.25$ and 
${\textit{Re}}=180$ at different phases: (a–d) in steps of 
$T/4$ with an exact Womersley profile imposed at the inlet. The flow is described by mode C.
Comparison of velocity profiles at inlet (plane of expansion) between exact Womersley and modulated Poiseuille for 
$U_r=1.25, a=0.25, Re=160$. Results are shown for (a) 
$t/T=0$, (b) 
$t/T=0.25$, (c) 
$t/T=0.50$, (d) 
$t/T=0.75$, (e) time averaged.
Streamlines for 
${{U_{r}}}=1.25$ at 
${{a}}=0.25$ and 
${\textit{Re}}=160$ at different phases: (a–d) in steps of 
$T/4$ with an exact Womersley profile imposed at the inlet. The flow is described by mode D.
For mode B, as discussed, there is splitting and diffusion of recirculation regions at a particular axial location. In the asymmetric regime, however (figure 22), the necking occurs at two axially distinct locations, dissimilar to the stable case. This results in waviness of the shear layer at the mid-section of the channel. The velocity perturbations (
$\tilde {u}$, 
$\tilde {v}$) show regions of prominent action (
$x=16-24$, figure 22) that results in the wavy nature of the shear layer. Here, we can see a number of perturbation recirculation eddies of the same sign near both the walls, resulting in asymmetric undulations in the shear layers. The Floquet multiplier for this case is 1.12 that is higher than 1, representing an unstable symmetric flow at this Re.
In the case of mode C, as discussed, the splitting of the recirculation region is present, but the distal eddy stays and diffuses in subsequent flow cycles. This results in the flow being more wavy in nature, downstream of the eddies (
$x=12-26$, figure 23). The action of the prominence of the velocity perturbations (
$\tilde {u}, \tilde {v}$) is downstream (
$x=12-26$) rather than near the expansion (
$x=2-8$). This may signify the occurrence of a convective instability. The perturbation vorticity field also shows the convective nature.
In mode D we observed an anti-phase response of the eddies, 
${{X_{r}^{*}}}$, with the incoming pulse, 
$u_m(t)$. The action of prominence is near the expansion (
$x=2-10$, figure 24). This may signify the existence of absolute instability. We further see counterclockwise long perturbation eddies over both the bottom and top walls.
In general, we can observe that larger perturbation eddies result in a lower critical Reynolds number (mode A and D), whereas shorter span multiple perturbation eddies result in asymmetric necking or splitting of recirculation regions with an increased critical Reynolds number.
Table 7 shows the Floquet multipliers obtained in each of the flow modes (A, B, C and D) along with 
${{Re_{c}}}$ from the Floquet study corresponding to 
${{Re_{c}}}$ obtained from our simulations. The overall trend matches well when variations are considered across different modes.
5. Conclusions
We have performed a detailed analysis of the effect of inflow pulsation on flow characteristics in a sudden expansion channel with a wide range of 
${{a}}$ and 
${{U_{r}}}$. In this regard, we have identified four different modes, namely, synchronized growth of the recirculation region (mode A), necking and diffusion of the recirculation region (mode B), splitting and convection of the recirculation region (mode C) and inverse growth of the recirculation region (mode D) corresponding to a very high to low time period of the pulsation cycle. For symmetry breaking, we observe the highest stability points for modes B and C and lower stability points for mode D. Symmetry breaking of mode A is lower than mode B and is similar to that of steady sudden expansion flows for 
${{a}}=0.25$. We have reported the mode map and found that at lower 
${{U_{r}}}$, mode D is observed irrespective of the amplitude. However, at moderate or higher 
${{U_{r}}}$, mode C and mode B are respectively promoted by increasing 
${{a}}$. At lower 
${{U_{r}}}$ (or time period), 
${{Re_{c}}}$ is smaller than the steady case, whereas it is more at higher 
${{U_{r}}}$. The departure from the steady value increases with 
${{a}}$. For a higher expansion ratio, 
${{Re_{c}}}$ is usually lower than the smaller 
${{ER}}$ case, as observed for steady flow. However, the trend is similar to the lower 
${{ER}}$ case. The convection and diffusion rates of vorticity are correlated with the vorticity dynamics of individual modes. Floquet analysis is carried out to estimate the stability of symmetric base flows and to explore the dynamics of infinitesimal perturbations in the symmetry breaking of the configuration. The obtained multipliers in each of the flow modes validate the 
${{Re_{c}}}$ from our simulation results.
Acknowledgements
The authors would like to thank Mr R.R. Prakash for his support in the current work. The authors acknowledge the use of the PARAMShakti supercomputing facility of IIT Kharagpur established under the National Supercomputing Mission (NSM), Government of India, and supported by the Center for Development of Advanced Computing (CDAC), Pune, India.
Funding
The present research is funded through SERB Project: CRG/2019/00265.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in Mendeley Data at https://doi.org/10.17632/ptzxkwyt2j.1.
Appendix A. Comparison with a long inlet channel upstream of expansion
In this section we compare our results with that of the inflow boundary condition (modulated Poiseuille, (2.1)) prescribed at the inlet of a long channel well upstream of the expansion section. We consider a channel of length 
$20$ and unity height upstream of the expansion plane. The grid size employed is similar to that obtained from the convergence tests performed in § 4.1. Firstly, we examine the critical Reynolds number for steady Poiseuille inflow. Table 8 reports the 
$Re_c$ for a long inlet channel upstream of expansion for 
$ER=2$, and is compared with Drikakis (Reference Drikakis1997). A good match is obtained. However, it is lower than the case of directly prescribing inflow boundary conditions at the expansion section (
${{Re_{c}}}=169$, table 4). Hence, the long inlet channel offers lower stability when compared with directly prescribing inflow at the expansion for the same 
${{ER}}$. This is also reflected in the pulsatile flow cases where the flow features are unstable and asymmetric in nature for the same parameters considered in § 4.2 for individual modes. Please note that the modes obtained are the same, but we obtain unstable asymmetric flow in the case of a long inlet channel for the same flow parameters. Similar to that of steady inflow (where a long inlet channel offers lower 
${{Re_{c}}}$ than directly specifying inflow at expansion), the 
${{Re_{c}}}$ in pulsatile cases might also follow the same trend, offering unstable asymmetric fields for the same flow parameters. Flow features of the individual modes are shown in figures 25, 26, 27 and 28 for modes A, B, C and D, respectively.
Appendix B. Comparison of results with Womersley profile imposed at plane of expansion
In order to ensure the generality of our results, we compare the same with the exact Womersley profile as inlet conditions (at the plane of expansion) instead of the modulated Poiseuille profile (2.1). Considering a driving pressure gradient
\begin{equation} -\frac{1}{\rho} \frac{\partial p}{\partial x} = p_0 + K \exp(\iota \omega t), \end{equation}
where 
$\iota =\sqrt {-1}$, 
$p_0$ is the steady part of pressure and 
$K$ is the amplitude of pressure oscillation. The resulting velocity profile reads
$$\begin{gather} u=3/2(1-4y^2)+3K/k (A_n(\cos(\omega t)+\sin(\omega t)) + B_n(\cos(\omega t)\nonumber\\ -\sin(\omega t)) +\sin(\omega t) - \cos(\omega t)), \end{gather}$$
where 
$\omega =2 {\rm \pi}/ T$, 
$k=\omega /(4 {\textit{Re}})$, 
$A_n$ and 
$B_n$ are given by
\begin{gather} A_n=(\cos(2y\beta)\sin(\beta)\cosh(2y\beta)\sinh(\beta)-\sin(2y\beta)\cos(\beta)\sinh(2y\beta)\cosh(\beta))/\varDelta, \end{gather}
\begin{gather}B_n=(\cos(2y\beta)\cos(\beta)\cosh(2y\beta)\cosh(\beta)+\sin(2y\beta)\sin(\beta)\sinh(2y\beta)\sinh(\beta))/\varDelta, \end{gather}
\begin{gather}\varDelta=\cos(\beta)^2 \cosh(\beta)^2 + \sin(\beta)^2 \sinh(\beta)^2, \end{gather}
where 
$\beta =\sqrt {k/2}$.
Figure 29 shows the comparison of velocity profiles at the inlet between modulated Poiseuille (2.1) and exact Womersley (B2) at four different instants of time and at temporal mean for 
$U_r=160, a=0.25, Re=170$, mode A. Both profiles closely match except for at 
$t/T=0$ (figure 29a) and 
$t/T=0.5$ (figure 29c). However, the mean flow profiles closely match (figure 29e). The flow features are also similar (figure 30 as compared with figure 3).
For mode B, we observe a slight deviation in the phase-wise distribution of velocity profiles (Womersley and modulated Poiseuille, figure 31a–d), but the mean profiles are similar (figure 31e). The flow field shows an overall delayed response due to the phase lag of the sectional mean Womersley as compared with that of modulated Poiseuille (compare the time chart of the flow present at the right-hand side of figures 32 and 6). A similar necking phenomenon is seen with prescribing the inlet profile as Womersley instead of modulated Poiseuille.
Similar deviations in phase-wise profiles are seen in the case of mode C (figure 33a–d), with the overall mean profile being similar (figure 33e). The features of splitting and advection of the recirculation region (mode C) is prevalent with the Womersley inlet, albeit by a phase lag (figure 34 as compared with figure 8).
For mode D also, we obtain changes in the phase-wise velocity profiles (figure 35a–d), while the flow features are similar, albeit by a phase lag (figure 36 as compared with figure 10).
Appendix C. Grid convergence study parameters
The order of convergence, 
$P$, is given by
\begin{equation} P=\frac{\ln \dfrac{F_3-F_2}{F_2-F_1}}{\ln r}, \end{equation}
where, 
$F_3$, 
$F_2$, 
$F_1$ are the length of the recirculation region at fine, medium and coarse grids and 
$r$ is the ratio of grid spacing between two successive grids. The predicted value of the recirculation length via Richardson's extrapolation, 
$REV$, is given by
\begin{equation} REV = F_3 + \frac{F_1-F_2}{r^P-1}. \end{equation}
The grid convergence index 
$({rm GCI})$ is given as
\begin{equation} {\rm GCI}=\frac{|e|}{r^P-1}, \end{equation}
where 
$|e|$ is the relative error between two successive grids. Lastly, the asymptotic range of convergence 
$(C)$ is obtained as
\begin{equation} C=\frac{{\rm GCI}_{2,3}}{r^P {\rm GCI}_{1,2}}, \end{equation}
where 
$({\rm GCI}_{2,3})$ and 
$({\rm GCI}_{1,2})$ are based on the medium-to-coarse and fine-to-medium grids, respectively.
Open access
