Periodic fluctuations of streamwise vortices in inertia-dominated intersecting flows

Vortices are a common feature in geophysical flows and are important for numerous industrial processes and engineering applications.^1,2 In the proximity of obstacles and stagnation points, vortical flows arise, and often become unsteady, at relatively low Reynolds numbers [ratio between inertial and viscous forces, Re ∼ O (10–100)].³ 

The unsteady motion of vortices can induce vibrations that are potentially destructive for structures (e.g., obstacles) such as bridges, pipes, buildings, and cables in both marine and terrestrial environments.^3,4 Vortex-induced vibration (VIV) is one of the major challenges in the engineering of deep water risers for oil and gas production, and many efforts have been made to suppress its effects.⁵ In order to avoid vibrations of underwater pipes, for example, they are covered with small staggered structures that are designed to stabilize the flow by suppressing vortex formation.^6–8 Additionally, VIV is a potential source for renewable energy from fluid flows,^9,10 particularly from ocean currents at high Re ∼ O (10⁴).¹¹ 

In recent years, there has been a growing interest in the additional influence of conveying fluids on the stability of pipes and risers.^12,13 It was found that above a critical internal flow velocity (constant or pulsating), VIV could match the frequency of an external flow-induced vortex shedding such that they resonate together, resulting in a premature failure of the riser system.^14,15 Similar internal vortex flow dynamics are also seen in inertial microscale flows.¹⁶ Vortex formation and fluctuations are important in inertial microfluidic applications such as cell sorting and stretching,^17,18 drug delivery through vortex hydroporation,¹⁹ mixing enhancement,^20–23 heat transfer enhancement,²⁴ and micro-oscillators.²⁵ 

The precise control over internal and external vortical flows, which arise in the proximity of stagnation points, can, therefore, be useful for improving the stability of structures that are subjected to currents, including pipes and channels conveying fluids at various length scales.

A possible approach to tune VIV is by introducing physical obstacles in the flow field. The dynamics of the vortical flow field and its feedback on the obstacle will depend on the shape, flexibility, and number of obstacles in the flow field as well as on the Reynolds number. Downstream of a bluff obstacle, spanwise vortices appear upon a critical Reynolds number, Re_c, due to the existence of a pinned stagnation point in the wake (e.g., the von Kármán vortex street).^26,27

Stagnation points are also present in intersecting flows (such as T-²⁸ and X-junctions with various crossing angles and aspect ratios^21,29,30), without the presence of physical obstacles. Four-way intersections such as the 4-roll mill and the cross-slot geometry with two opposing inlets and two opposing outlets produce a “free” stagnation point, formed by the conveying fluid, that is not pinned at any fixed surface in the flow. In such cases, for Re > Re_c, the formation of axisymmetric streamwise vortices is observed, as illustrated in Fig. 1.^21,31,32

In axisymmetric vortices, the vorticity is concentrated by axial stretching and diffused by viscous stresses.³³ There are various known models that describe axisymmetric vortices, all of which are solutions to the Navier–Stokes equations. Such models include Lamb–Oseen, Burgers, Newman, Sullivan, Batchelor, and others.³⁴ Interestingly, all these vortex models have similar characteristics and can be grouped into a single family with common mathematical traits. These traits are captured by the sine-like velocity profile of the vortex cross section. More specifically, the Burgers and Lamb–Oseen models describe stagnation point vortices. The Burgers vortex is defined as a stagnation point stretched vortex for which the viscous diffusion is balanced by the concentration process induced by stretching. The Lamb–Oseen vortex is a non-stretched stagnation point vortex that is mainly affected by viscous diffusion. Nevertheless, the velocity profiles of the Burgers vortex and the Lamb–Oseen vortex are identical when normalized.³⁵ 

Experimental evidence of Burgers-like vortices (Kerr–Dold vortices) was found in extensional flows, induced by cross-slot geometries with high aspect ratios, α = 32 and α = 262 (α = d/w, where d and w are the depth and width of the geometry, respectively; see Fig. 1).^32,36 In these cases, the vortices emerge as a stack of counter-rotating streamwise vortex pairs upon Re_c ≈ 40 and Re_c ≈ 55 for the lower and higher aspect ratios, respectively. In these reports, the characteristic length scale used to determine the Reynolds number is the width of the channel, w; see Fig. 1.

For flows beyond a critical Reynolds number in the cross-slot geometry with a low aspect ratio (α ≲ 2.5), a single streamwise Burgers-like vortex forms at the intersecting region; see Fig. 1.^21,29 The transition from a symmetric flow field to an asymmetric flow field is a bifurcation and can result in either a clockwise (CW) or a counter-clockwise (CCW) rotating vortex. The transition occurs upon a value of Re_c that is α-dependent, with a tricritcal point at α = 0.55. For an aspect ratio α < 0.55, the transition will be supercritical, while for α > 0.55, the transition will be subcritical with a characteristic hysteresis. For increasing aspect ratio in the range 0.4 ≲ α ≲ 2.5, the onset Reynolds number reduces from Re_c ≈ 110 to Re_c ≈ 25.^21,29 Similarly, the characteristic length scale used to determine the Reynolds number is the width of the channel, w; see Fig. 1.

Identification of streamwise vortices (and vortices in general) in a flow field can be achieved in numerous ways; however, there is no universal acceptance of a single definition.^37–39 A few of the most common methods to identify vortices in planar velocity fields are based on the local velocity gradient tensor,⁴⁰ including the so-called λ_ci criterion⁴¹ that was recently used to identify vortices in a microfluidic mixing–separating device.⁴² Vortex identification by streamlines that exhibit circular or spiral patterns is another common tool to identify vortices in the flow field.⁴³ Although such a seemingly simple approach is not a robust method to identify vortices in turbulent flows, it can be applied at low to moderate Reynolds numbers, when the observation is done over the normal plane to the vortex core and the reference frame is steady.⁴⁴ 

Vortical flows may become unsteady when inertial forces are increased. For a certain range of Re, these flows may exhibit time-dependent periodic cycles such as the von Kármán vortex street that develops downstream of a bluff body [in two-dimensional (2D) and three-dimensional (3D) flows].³ This flow instability emerges upon a relatively low Re_c ≈ 10 (the characteristic length scale used to determine Re is the bluff body diameter), such as two attached counter-rotating, spanwise vortices flowing downstream of an obstacle. The vortices grow in size for Re > Re_c.³ When a second critical Reynolds number, Re_p ≈ 40, is reached, periodic flow patterns emerge as the vortices that are formed behind the obstacle are shed in an alternate manner, generating VIVs that are exerted on the obstacle.^3,45

Periodic streamwise vortex flows have also been predicted by numerical simulations at flow intersections and have recently been experimentally validated in intersections such as the T-shaped geometry. This geometry is well studied for its mixing efficiency^46,47 and rich flow regimes observed due to the generation of various vortical flows, namely, symmetric Dean-vortex flow, asymmetric vortex flow (e.g., engulfment flow), and unsteady flow regimes.⁴⁸ This periodic motion can be characterized by its frequency, f, or with a more general non-dimensional frequency, the Strouhal number, St. The typical onset of the unsteady periodic flow regime was found at Re ≃ 200, with St ≃ 0.1.^23,47–51

Periodic flows with similar dependence between St number and Re were also observed in a microfluidic T-channel with offset inlets.⁵² The periodically pulsating state of bubble-type vortex breakdown was reported at 120 ≤ Re ≤ 170, St ≃ 0.5, while a helically oscillating state of the spiral-type vortex breakdown was found at Re ≥ 160, St ≃ 1.75.

Similar to the flow regimes that were reported at cylinder wakes and at T-intersections, in the cross-slot geometry, it was also observed that the streamwise vortex would grow in size for Re > Re_c.^21,29,53 However, the transition to unsteady flow has only been reported in one recent study at a single value of the aspect ratio, α = 1.²² In this study, dye advection patterns were used to qualitatively analyze the flow fluctuations. Additionally, the same authors have recently reported that trapping of dye at the center of the vortex region occurs due to a bubble-type vortex breakdown.⁵⁴ Quantitative measurement, theoretical description, and numerical solution of the unsteady vortex flow regime remain to be explored.

Here, we use quantitative time-resolved flow velocimetry measurements, supported by 3D numerical simulations, to characterize steady vortical flow fields in two cross-slot geometries with different aspect ratios, α. We show that by changing α, it is possible to precisely modulate certain features of the steady-state vortices, which are generated for flows with Re > Re_c. We demonstrate control over features including the number of vortices in the flow field, their relative rotation direction, and vortex core structure. We show that the unique steady-state features determine the onset parameters of time-dependent flows for Re > Re_p and the nature of the fluctuations that arise. Additionally, we find that the relation between Re and St numbers, which corresponds to the different modes of the fluctuations, can be described with the same quadratic relation that is found to occur for external flows around circular cylinders.

Two microfluidic cross-slot devices with approximately inverse values of the aspect ratio are selected for this study. The devices are fabricated by selective laser-induced etching in fused silica using a “LightFab” 3D printer (LightFab GmbH, Germany).^55,56 These specific devices are selected due to their positions at the extremes of the Re–α phase diagram [see supplementary material file (SI), Fig. S1]. These devices enable comparison between the different types of vortical flow fields that may form at the stagnation point of a four-way intersection.

The flow field is directly controlled by the Reynolds number, describing the ratio of inertial to viscous forces in the flow,

where U = Q/wd is the average flow velocity in each inlet and outlet of the cross-slot and ν is the kinematic viscosity of the fluid (Q is the imposed volumetric flow rate in each inlet and outlet). All of the experiments are done with deionized water (ν = 8.9 × 10⁻⁷ m² s⁻¹) at 25 °C.

The first geometry is designed with a depth d = 1500 μm and a width w = 620 μm, forming a relatively high aspect ratio of α = 2.4 [Fig. 2(a)]. In this channel, a subcritical transition, with hysteresis, from a symmetric to an asymmetric flow occurs at Re_c ≈ 26. Further details about this transition can be found in our previous publications.^21,29 The ratio between the vortex radius and the channel depth is 2R/d ≈ 0.5 for 32 ≤ Re ≤ 62, where R is the radius of the vortex, defined where v_y = v_y,max, at the coordinate (x/w, y/w, z/w) = (0, 0, −0.6), as shown in Figs. 2(b) and 2(c). As evident from experiments and numerical simulations, the vortex that is formed for α = 2.4 occupies a larger space than the downstream width, w, of the outlet channel [Figs. 2(d)–2(f) (Multimedia view)]. The borders of this vortex may intrude into the inlet channels (along the y-direction), and therefore, w is not a restricting factor for the vortex radius at the x = 0 plane.^21,29

The second geometry is designed with depth d = 670 μm and width w = 1490 μm, forming a relatively low aspect ratio, α = d/w = 0.45 [Fig. 3(a)]. A supercritical transition from the symmetric to the asymmetric flow with a central vortex occurs at Re_c = 107.²⁹ For a steady flow field and Reynolds number range 400 ≤ Re ≤ 446, we find that the ratio 2R/d ≈ 0.44. Here, R is defined where v_y = v_y,max, at the coordinate (x/w, y/w, z/w) = (0, 0, −0.1); see Figs. 3(b) and 3(c). The central stretched streamwise vortex that is formed for the case of α = 0.45 has minimal interaction with the borders of the outlet channel as indicated by experiments and numerical simulations [see Figs. 3(d)–3(f) (Multimedia view)].

The inlet lengths of the devices are set to 16 mm, providing a high ratio (>10) between the inlet length and the larger channel dimension and ensuring a fully developed flow before the fluid reaches the intersection region of the cross-slot geometry. The maximal development length of the flow for the highest Re studied (α = 0.45 and Re = 543) is $lh,laminar=0.05DhReDh=15$ mm, where D_h = 2wd/w + d = 0.9 mm and $ReDh=UDh/ν=337$⁠. The outlets of the channel are designed to be as long as possible (≃5 mm), while allowing imaging to be performed at the yz center plane (where x = 0).

The flow is driven using four high precision neMESYS syringe pumps (Cetoni GmbH, Germany), each fitted with a Hamilton Gastight glass syringe with a volume of 25 ml. Two of the pumps drive the fluid into the two opposing inlets, while the other two pumps simultaneously withdraw the fluid from the two opposing outlets (all syringe pumps are set to an equal volumetric flow rate, Q). The pumps are operated at a minimum of 10× (and typically >50×) the manufacturer’s specified lowest “pulsation-free” dosing rate. The system’s compliance is reduced to a minimum by using rigid poly(tetrafluoroethylene) tubing connections between the syringes and the microfluidic devices. The tubing is kept as short as possible, and all air bubbles are carefully purged from the system prior to all experiments.

Quantitative measurements of the flow field are done with a micro-particle image velocimetry (μ-PIV) system (TSI, Inc., MN).^29,53 Fluorescent tracer particles are seeded in the fluids for flow visualization (PS-FluoRed, MicroParticles GmbH, Germany). The diameter of the particles is d_p = 5.0 μm with excitation and emission wavelengths of 530 nm and 607 nm, respectively. The microfluidic devices are mounted on the stage of an inverted microscope (Nikon Eclipse Ti), and a Nikon Plan Fluor, 4×, NA = 0.13 numerical aperture lens is used to bring the central x = 0 plane into focus. With this combination of particle size and objective lens, the measurement depth over which particles contribute to the determination of the velocity field is δx_m ≈ 210 μm (δx_m = 0.33w for α = 2.4 and δx_m = 0.14w for α = 0.45).⁵⁷ Pairs of particle images are captured using a high-speed 1280 × 800 pixel CMOS camera (Phantom MIRO) operated in frame-straddling mode and synchronized with a dual-pulsed Nd:YLF laser light source with a wavelength of 527 nm (Terra PIV, Continuum, Inc., CA). The time separation between laser pulses is set such that the average particle displacement between the two images in each pair is ∼4 pixels. Cross correlation of each image pair is done with a standard μ-PIV algorithm (TSI, Inc.) that produces velocity vectors with components v_y and v_z in the y- and z-directions, respectively. For each frame, the vectors are spaced on a square grid of 32 μm × 32 μm in y and z. These vectors are then used to calculate and plot vorticity (ω_x) contours across the x = 0 plane using a Matlab code. Additional specifications about the μ-PIV system, the measurement, and processing methods can be found in our previous publications.^29,53

Using the syringe pumps to control the flow rates, a steady vortical flow is induced in the cross-slot by setting the Reynolds number to Re_c < Re < Re_p, where Re_p is the critical Reynolds number for the onset of time-dependent flow. Time-dependent flows are then induced by setting Re > Re_p. In order to avoid measuring transients, Re is kept constant for 5 s prior to the beginning of data acquisition. Note that this is a relatively long time in terms of convective time scale, t_c = w/U = 0.013 s for the lowest imposed flow rate studied here (α = 2.4 and Re = 32). Image pairs are captured at a high frequency of 500 Hz in order to resolve the details of time-dependent periodic cycles in the flow (apart from a measurement taken for α = 0.45, Re = 543, where the acquisition rate was 1000 Hz). For each Re studied, 100 velocity vector fields are produced in order to (1) confirm that the flow is steady for Re < Re_p and (2) capture several cycles of the periodic fluctuations for Re > Re_p.

The computational fluid dynamics (CFD) simulations performed in order to accompany the experimental findings consider a laminar, incompressible, and isothermal Newtonian fluid flow. The continuity and the momentum equations are solved numerically considering the finite-volume method,⁵⁸ 

where u is the velocity vector, ρ is the fluid density, p is the pressure, t is the time, and τ corresponds to the stress tensor. The discretized partial differential equations from Eqs. (2) and (3) are solved using an in-house implicit finite volume CFD solver, developed for collocated meshes^59,60 considering a generalized coordinate system. The pressure and velocity fields are coupled using the SIMPLEC algorithm for collocated meshes by employing the Rhie and Chow interpolation technique.⁶¹ In Fig. S2 of the supplementary material, the independence of the numerical solution to the numerical mesh is demonstrated. The convection terms in the momentum equation are approximated using the CUBISTA high resolution scheme,⁶⁰ the diffusive terms are discretized considering a central-difference scheme, while the transient term in the momentum equation is evaluated using a first-order implicit Euler scheme. The length of the two inlet channels is set to be 10w. To ensure that the flow is fully developed before approaching at the intersection of the channels, the fully developed velocity profiles are applied at the inlets of the domain. The profiles are calculated based on the value of the aspect ratio of the rectangular cross section,⁶² 

where a and b are the half-width and half-depth of the channel cross section, respectively. At the outlet boundaries, zero streamwise gradients are applied,⁵⁹ while the lengths of the outlet channels are set to be 50w, which is long enough to ensure that potential numerical artifacts from the outlet boundaries do not influence the flow upstream. Finally, the numerical results presented for both geometries are those obtained with the most refined meshes considered. More details can be found in Table S1 in the supplementary material.

For the aspect ratio α = 2.4, the transition to a steady vortical flow is subcritical, and the central streamwise vortex grows and intensifies at relatively low Re [see Fig. 2 (Multimedia view)]. Additionally, as mentioned in Sec. II A 1, the width of the channel (w) is not a confining factor since the vortex can grow beyond y/w = ±0.5. At an intersection with α = 2.4, a single steady vortex will form at Re_c = 26.5, as previously reported.^21,29 For this aspect ratio, the steady central streamwise vortex persists over a range of 26.5 < Re < 69. As described in the Introduction, for flows that are not turbulent (when the observation is done over the normal plane to the vortex core and the reference frame is steady), a simple streamline identification scheme can be applied to locate the vortices in the flow field. A contour plot of the normalized velocity, v_y/U, from the experimental results, is constructed for Re = 50 [see Fig. 2(b)]. The velocity profile along y = 0 is plotted as a function of z/R. The vortex radius R is determined at the coordinate at which v_y/U = v_y,max/U, as shown in Fig. 2(c).

Contour plots of the normalized vorticity are constructed for Re = 37 and Re = 50 from numerical simulations [see Figs. 2(d) and 2(e)] and experimental results [see Fig. 2(f) (Multimedia view)]. The flow field is dominated by a CCW vortex that is located at the center of the x = 0 plane. Numerical simulations reveal the 3D structure of the flow field at Re = 37 and Re = 50, showing the stretching of the vortex, by the streamwise velocity gradient, from the stagnation point located at (x, y, z) = (0, 0, 0) in opposing directions toward the outlets of the channel [the x-direction, Fig. 2(d)]. Please note that the vortex is stretched in both outlet directions, even though the perspective view slightly hides the stretching in the negative x-direction. When a non-confined vortex is stretched, its cross sectional area is reduced; hence, the vorticity intensifies along the vortex to preserve the angular momentum. The stretching is balanced by viscous diffusion, which causes the dissipation of vorticity away from its core.³⁶ In our setup, the geometric confinement balances the stretching process leading to a steady vortical flow field for 26.5 < Re < 69.

The 2D view of the yz plane from the numerical simulations at the Reynolds numbers shown in Fig. 2(e) indicates the change in the vortex core structure, in agreement with the experiments shown in Fig. 2(f) (Multimedia view). Vorticity patches with CW rotation appear where z/w ≳ 0.9 and z/w ≲−0.9. From the numerical simulations, we can see that the vorticity patches seen for Re = 37 develop to small CW vortices at Re = 50 [see Fig. 2(e)]. Figure 2(g) illustrates the steady flow vorticity profile along the z-direction, extracted from experiments and numerical simulations for Reynolds numbers 32 ≤ Re ≤ 62. Both experiments and simulations show a clear difference in the generated response when Re is increased. For this range of Reynolds numbers, we can see a gradual development of the vorticity profile, from a profile with a single peak (for Re = 32) to a profile with a central dip at z = 0 that increases with increasing Re. The development of a ring of high vorticity surrounding a lower intensity vorticity region most likely occurs due to a bubble-type vortex breakdown^63,64 [see Fig. 2(d)]. A vortex breakdown of this type was recently reported to occur in a cross-slot geometry with α = 1 at Re = 190.⁵⁴ Vortex breakdown was also reported in dividing T-junction flows in microfluidic channels.^52,65 Further analyses of these vortex profiles are presented in Sec. III B. We note that there are no existing vortex models to describe the vortical flow field that is observed here.

For this relatively low aspect ratio, α = 0.45 [see Fig. 3 (Multimedia view)], the transition from a symmetric to an asymmetric steady vortical flow is supercritical with a relatively slow transition dynamics occurring at a critical Reynolds number, Re_c = 107. A single steady vortex will persist and remain steady for a range of 107 < Re < 473.

A contour plot of the normalized velocity, v_y/U, from the experimental results, is constructed for Re = 418 [see Fig. 3(b)]. The velocity profile along y = 0 is plotted as a function of z/R. The vortex radius R is determined at the coordinate at which v_y/U = v_y,max/U, as shown in Fig. 3(c).

A contour plot of the non-dimensional vorticity, ω_xw/U, is constructed from numerical simulations and experimental data for a steady vortex flow at Re = 418 [Figs. 3(d)–3(f) (Multimedia view)]. Numerical simulations reveal the 3D structure of the flow field at Re = 418, showing the stretching of the vortex in two opposing directions toward the outlets of the channel [x-direction, Fig. 3(d)]. A 2D view of the yz plane from the numerical simulations at Re = 418 [Fig. 3(e)] shows good agreement with the experimental results [Fig. 3(f) (Multimedia view)].

Additionally, in Fig. S3 of the supplementary material, the independence of the numerical solution to the numerical mesh is demonstrated.

In Figs. 3(e) and 3(f), it is seen that the flow field is dominated by a central CCW rotating vortex and two smaller co-rotating vortices at the top-left and bottom-right corners. Additionally, there are two CW rotating vortices located between the central vortex and the side vortices, in the top-left and bottom-right corners of the plane.

The normalized vorticity, ω_x/ω_x,max, obtained from the experiments and from the numerical simulations for Re = 418 is plotted as a function of z/R (x = y = 0) in Fig. 3(g). A Gaussian-like profile of the vorticity is seen at the yz plane along the z-direction for the range −1.5 ≲ z/R ≲ 1.5 (y = 0). The peak of the vorticity is at the location of the stagnation point. Similar vortices, with a single peak in the vorticity distribution and Gaussian-like profiles, that do not have counter-rotating vortices with similar intensity surrounding them are commonly seen in stratified fluids and are described as “monopolar” vortices.^66–68 These vortices are confined in both streamwise and spanwise directions.

Although a simple Gaussian model can describe the vorticity profile for the vortices formed for α = 0.45 [Fig. 3(g)], it will not be able to capture the complexity of the core structure that is seen for α = 2.4 [Fig. 2(g)]. In Subsection III B, we will present a simple modification to an existing model that describes pancake-like vortices,⁶⁶ and we will test its ability to describe the different vortex profiles seen here.

Confinement of vortices, within a cross-slot channel of the low aspect ratio, α, results in the formation of a single, isolated, and steady vortex while maintaining the streamwise component of the vortex, allowing it to stretch downstream (see Figs. 1–3). The description of such a vortex fits well with the Burgers vortex model; however, a confinement factor is missing from the Burgers description.⁶⁹ More specifically, the Burgers model (as well as Lamb–Oseen, Rankine, and others) does not account for the sharp decrease in the velocity profile that occurs in the case of confined streamwise vortices.³⁴ In these models, the vorticity profile is dominated by the balance between viscous diffusion and vortex stretching. However, in our experimental setup, the no-slip condition at the walls of the outlet channel imposes a forced velocity drop, as seen in the velocity profile v_y along the z-direction [Figs. 2(c) and 3(c)]. Therefore, viscous diffusion does not have a significant influence on the flow field. Additionally, these models always describe a peak in the vorticity distribution at the center of the vortex, which is not always the case in the current flow, as we show in Fig. 2(e). We have found no existing model that describes a single confined Burger-like vortex.

A model for confined vortices was previously reported for flows that are spatially confined within a stratified fluid layer.⁷⁰ These non-stretched vortices are confined in both streamwise and spanwise directions and resemble a pancake shape. Pancake-like vortices have a forced drop velocity profile, v_y along z (or v_z along y), due to their confinement. Consequently, viscous diffusion is not the main factor influencing the velocity profile of these vortices.⁶⁶ It is important to note that the forced velocity drop not only affects the far edges of the vortex but also the structure and stability of the vortex core.⁷¹ 

Here, we modify the pancake-like vortex model, suggested by Flór and van Heijst,⁶⁶ to describe the vorticity profiles of the vortices in the cross-slot geometry that are confined only in the spanwise direction (and transverse directions once in the outlets) and are stretched in the streamwise direction. There are several ways to define the vortex diameter, but in order to apply the pancake-like vortex model, we use the same definition as in Flór and van Heijst.⁶⁶ The vortex radius, R, is defined at the point where v_y = v_y,max along z at x = y = 0. To apply the pancake model for the vorticity distribution, ω_x along z is re-scaled with ω_x,max.

For α = 0.45, the peak of the vorticity distribution is at the center of the vortex [at z/R = 0, Fig. 3(e)]. However, for α = 2.4 (37 < Re < 50), the peak of the vorticity is not located at the center of the vortex [Fig. 2(e)]. Such complex vorticity profiles are not common, and we could only find scarce evidence in numerical studies of fluid drainage systems⁷² and for tangential wind profiles of vortices.⁷³ In order to capture the complexity of the vorticity profiles within the vortex core, the pancake-like vortex model⁶⁶ is adjusted as presented in Eq. (5),

where b₁, b₂, c₁, c₂, R₁, and R₂ are free parameters, and their best-fit values by the model above are specified in Table I. The parameters b₁ and b₂ are related to the orientation and intensity of the vorticity. For α = 0.45, where the vorticity has a maximum at the center, we get b₂ ≈ 0, canceling the second term of the equation. The radii of the vortex and the vortex core are depicted by R₁ and R₂, respectively.

The parameters c₁ and c₂ account for the steepness in the vorticity profile. Values of c₁ ∼ 2 indicate that a pancake-like vortex is close to the critical parameters in which it will become unsteady.⁶⁶ 

Figure 4 shows model fittings with Eq. (5) to numerical simulations and experimental results for α = 0.45 at Re = 418 [Figs. 4(a) and 4(b)] and α = 2.4 at Re = 37 [Figs. 4(c) and 4(d)] and Re = 50 [Figs. 4(e) and 4(f)]. In all these cases, the flow field is within the steady vortex flow regime. From the plots shown in Fig. 4, we can see that Eq. (5) describes well the monopolar central vortex at α = 0.45. With our added second term, the model also describes well the double peak vorticity profiles observed for α = 2.4. The model does not apply to the areas of counter-rotating vorticity close to the bounding walls, where ω_x/ω_x,max < 0.

Comparing the fitting parameters that are presented in Table I, we see the close agreement between experimental and numerical results. The parameters b₁ and b₂ show close agreement in all cases. The radius of the vortex, R₁, and the radius of the core, R₂, show slight shifts between experiments and numerical simulations. While the vortex radius R₁ remains fairly constant as Re is increased (for α = 2.4), we find that the radius of the core, R₂, grows with increasing Re. The increase in the radius of the vortex core occurs due to increasing inertia (increasing Reynolds number), which leads to a pressure rise along the center of the vortex core, in the streamwise direction, and results in an expansion of the vortex core.⁷⁴ The steepness of the vorticity profile of the vortex and the core slightly increases when Re is increased, indicating that the vortex is close to the critical value in which it will become unsteady (c₁ ∼ 2). For the measurements conducted here, the vortices remain steady for values of c₁ > 2. This can be explained by the increase in the axial stretching rate (induced by an increase in Re) that was found to have a stabilizing effect on non-confined Burgers vortices.³⁵ Axial stretching does not occur in pancake-like vortices, and therefore, they become unstable at lower values of c₁.

The parameter c₂ increases when the Reynolds number is increased from Re = 37 to Re = 50, indicating that the dip of the vorticity in the vortex core grows deeper. We note that the numerical simulations result in a deeper local minimum at the center of the vortex core, which can be explained by the experimental measurement depth (δx_m ≈ 0.33w) and the higher spatial accuracy of the resolved velocity fields of the numerical measurements.

Time-resolved PIV experimental measurements and numerical simulations reveal that time-periodic flows emerge in the cross-slot geometries when imposing Re > Re_p. To characterize these flows, we define a dimensionless time t^*,

For the case of α = 2.4, the flow field becomes unsteady at Re > 69. For α = 2.4, the velocity component v_y is extracted from each measurement from the coordinate (x/w, y/w, z/w) (0, 0, −0.6), at which the fluctuations in the flow are the strongest and the detected signal is particularly clear. At each measurement, Re is kept constant for a range of 50 ≤ Re ≤ 249.

A time series of instantaneous contour plots of ω_xw/U from experiments and numerical simulations for Re = 87 is shown in Figs. 5(a) and 5(b) (Multimedia view), respectively. Here too, the experiments are conducted at 500 Hz, with a time step of 0.002 s between two images; a full cycle of the fluctuation at Re = 87 is completed within 0.022 s. The period time found in both experiments and numerical simulations is t^* = 4.0.

The nature of the fluctuation in both experiments and numerical simulations can be described as follows: from t^* = 0 to t^* = 1.6, the central vortex maintains a large core with a distinct ring of intense vorticity around the core, the core of the vortex slightly grows, and at t^* ≈ 2.0, the core structure begins to deform [Figs. 4(a) and 4(b)]. At t^* = 2.4 and t^* = 2.8, the vortex core shrinks to a smaller, concentrated structure with an intense vorticity at its center. At t^* ≈ 3.2, the core of the vortex is reformed, and the vorticity at the center of the vortex becomes less intense than the surrounding ring. At t^* = 3.6 and t^* = 4.0, the vortex core begins to grow. From the numerical simulations, we can also see that the vorticity patches located at z/w > 0.9 and z/w < −0.9 develop to CW rotating vortices and migrate toward z/w = 0, where they interact with the central CCW vortex as it collapses [see Fig. 5(b) (Multimedia view)]. The changes in the vorticity distribution of the central CCW vortex, as the core collapses and reform, lead to changes in the local shear-stress and to a response from the surrounding CW rotating side vortices that follow the central vortex movement.⁷⁵ 

In Fig. 5(c), the vorticity distribution along the center line (z/w) is plotted. These are taken from the experimental contour plots shown in Fig. 5(a) (Multimedia view). From t^* = 0 to t^* = 2.0, a double peak is seen, indicating an intense vorticity ring surrounding the core of the vortex. At t^* = 2.4, the double peak is replaced by a single peak, indicating that the vortex core has lost its structure. At t^* = 2.8, a single peak is seen in the vorticity distribution, and at t^* = 3.2, the double peak reappears, indicating that the vortex core and the surrounding ring of the intense vorticity region has been re-formed. The fluctuations of the vortex core are attributed to a bubble-type vortex breakdown that appears in both outlet arms. A similar flow was recently reported to occur in a cross-slot geometry with an aspect ratio α = 1 at Re = 190.⁵⁴ In Fig. 2(d), there is a clear indication that a bubble-type vortex breakdown occurs. We also observe that the vortex breakdown bubble expands further downstream when the Reynolds number is increased from Re = 37 to Re = 50.

For each Reynolds number studied, the velocity component v_y|_0,0,−0.6 is normalized by subtracting the time-averaged velocity $v¯y|0,0,−0.6$⁠, according to Eq. (7), such that the data fluctuate around v_y = 0,

A comparison of the periodic signal found in experiments and numerical simulations is presented in Fig. 5(d). The normalized re-scaled velocity, v_n, is plotted as a function of time, t^*, for α = 2.4 and Re = 87. The period time of the fluctuation and the magnitude of v_n show close quantitative agreement between experiments and simulations.

In order to extract the characteristic frequencies of the fluctuations, the power spectral density (PSD) is found by applying a fast Fourier transform (FFT) to the normalized velocity time-resolved series v_n [determined according to Eq. (7)], plotted in Fig. 5(d). We define the non-dimensional frequency, the Strouhal number, St, as

where f is the frequency of the fluctuation. We find that using PSD reduces the noise in the measurement and produces a clear signal from which we could extract the typical frequency and calculate the Strouhal number. In Fig. 5(e), the PSD is plotted as a function of St.

FFT analysis is conducted on the data presented in Fig. 5(d). In Fig. 5(e), the PSD is plotted as a function of St, showing quantitative agreement between the experiments and simulations in the fundamental frequency, f₀, of the fluctuations as well as the harmonics.

Additional experiments are performed for a range of 32 ≤ Re ≤ 249; see supplementary material, Fig. S4. According to these measurements, f₀ increases as Re increases; however, St behaves in a non-linear manner. These findings will be further examined in the discussion in Sec. IV.

A time series of the velocity component v_y for the case α = 0.45 is extracted from each measurement from the coordinate (x/w, y/w, z/w) (0, 0, −0.1), where the fluctuations in the flow are the strongest and the detected signal is particularly clear. At each measurement, Re is kept constant in the range 445 ≤ Re ≤ 543. Due to limitations of the current experimental setup, we could not investigate higher Re, as our system is limited by the volume of the syringe pumps and the capabilities (i.e., maximum resolvable velocity) of the μ-PIV system.

For the case of α = 0.45, when Re is increased to Re ≥ 473, the flow field loses steadiness and periodic fluctuations emerge. A high temporal resolution μ-PIV measurement of the flow field (at 500 Hz, time step between images is 0.002 s) reveals that the fluctuations have a period time of 0.016 s (t^* = 2.8), for 473 ⩽ Re ⩽ 543.

A time series of instantaneous contour plots of ω_xw/U are constructed from experimental and numerical data for Re = 488 and shown in Figs. 6(a) and 6(b) (Multimedia view), respectively. From the experimental results [Fig. 6(a) (Multimedia view)], the nature of the fluctuation can be described as follows: at t^* = 0, the flow field contains a large CCW rotating vortex, with its axis located at y = z = 0. Additionally, two smaller, less intense CCW rotating vortices can be seen in the upper left (x/w, y/w, z/w) ≈ (0, −0.35, 0.15) and lower right (x/w, y/w, z/w) ≈ (0, 0.4, −0.15) corners of the yz plane. From t^* = 0 to t^* = 2.0, the upper left corner vortex gradually intensifies as it migrates toward the center of the yz plane; simultaneously, the central vortex loses intensity and migrates toward the bottom right corner of the yz plane. At t^* = 2.0, the two vortices begin to merge as they share thin lines of concentrated vorticity (i.e., vortex filaments⁷⁶) and rapidly move toward each other. At t^* = 2.8, the vortices have completed the merging process as they form a single vortex at y = z = 0. A new small vortex forms at the upper left corner of the plane [Fig. 6(a) (Multimedia view)].

Snapshots of the periodic cycle found by the numerical results [Fig. 6(b) (Multimedia view)] show a more complex behavior, although some similarities can be observed when compared with the experimental results. The nature of the fluctuation can be described as follows: from t^* = 0 to t^* = 3.15, the top-left side vortex intensifies and simultaneously migrates with the central vortex toward the bottom-right corner, resulting in vortex merging at t^* = 3.85. This part of the fluctuation is similar to the full fluctuation cycle found in the experimental results [Fig. 6(a) (Multimedia view)]. Following the completion of the first vortex merging process from t^* = 5.25 to t^* = 7.00, a second merging process is seen; however, now it is the bottom-right side vortex that intensifies and migrates simultaneously with the central vortex toward the top-left corner of the plane, leading to vortex merging at t^* = 7.35.

In the experimental results, it is always the top-left vortex that intensifies and merges with the central vortex, while in the numerical simulations, both the side vortices intensify and merge with the central vortex in an alternate manner. This divergence may be related to inherent small asymmetries present in the experimental microfluidic channel that favor the intensification of the top-left side vortex. By contrast, in the numerical simulations, which are performed in a perfectly symmetric device, there is no favorable side-vortex that intensifies.

In Fig. 6(c), the vorticity distribution ω_xw/U along y = 0 from the time series shown in Fig. 6(a) (Multimedia view) is plotted as a function of z/w. The cycle starts at t^* = 0 with a peak at z/w = 0, where ω_xw/U > 50, and the vortex is at the center of the yz plane. The peak intensity is gradually reduced to ω_xw/U < 10, where the central vortex is at the farthest point from the center at t^* = 1.6. At t^* = 2.0, ω_xw/U along z/w rapidly increases as the vortices merge at the center of the yz plane.

The signal of the velocity v_y|_(0,0,−0.1) found in experiments and numerical simulations is compared for the case where α = 0.45 and Re = 488 [Fig. 6(d)]. In Fig. 6(e), the PSD is plotted as a function of St. It is visible that the signals of experiments and simulations occasionally overlap. However, the velocity time series extracted from the numerical simulations does not result in an ordered periodic signal with a single dominant frequency.

Interestingly, we observe in Fig. 6(e) that the fundamental frequency f₀ obtained from the experimental results coincides with the fourth harmonic found from the numerical simulation. The overlap between these peaks can be attributed to the merging of the top-left vortex with the central vortex, which is observed in the results from both experiments and simulations.

To exclude the case of a numerical artifact due to the employed time step, we provide in Fig. S5 of the supplementary material a comparison between different time steps, where it can be seen that good agreement exists. According to these plots, we see that the numerical simulations result in a semi-periodic fluctuation in all time steps that were studied. Additionally, three other cases are investigated numerically to find the effects of the boundary conditions (short outlet length, uniform inlet velocity profile, and the combination of both) on the flow field; these are shown in Fig. S6. It is found that the fluctuation remains semi-periodic in the cases where the outlet is short (4.5w) and the case in which the outlet is short and the inlet flow profile is uniform. For the third case, in which only a uniform inlet flow profile is imposed (outlet length is 50w), the extracted signal shows an ordered periodic behavior. Thus, the precise dynamics are sensitive to the boundary conditions used in the simulations. Aside from slight geometric biases which may be present in the experiments, as already discussed, the geometric setup required of the experimental domain to obtain μ-PIV measurements necessarily required finite inlet and outlet channels, and these may be the cause of these subtle differences observed with simulations. Future studies where these lengths are systematically varied should help shed light on this topic, but are beyond the scope of the current study.

Additional experimental measurements are performed for 473 ≤ Re ≤ 543, and the plots of PSD as a function of frequency, f, can be found in the supplementary material file, Fig. S7. According to these measurements, f₀ remains stable in the range 473 ≤ Re ≤ 543; hence, St is decreased (as increasing Re is here achieved by increasing the flow rate). These findings will be further examined in Sec. IV.

The natures of the fluctuations that are observed in the two different cross-slot channels studied here (α = 2.4 and α = 0.45) are fundamentally different due to the initial flow field that is dictated by the aspect ratio. For α = 2.4 (supercritical transition from the symmetric to the asymmetric flow), the steady flow field, prior to the onset of the periodic flow, is dominated by a single CCW rotating vortex with a core structure that resembles an “eye” (the center of the core has a dip in the vorticity). Additionally, above and below the central vortex, there are two CW rotating vorticity patches that develop to vortices during the fluctuation. Although the fluctuation does not exhibit vortex merging or splitting, the central vortex core structure collapses and reforms, while the smaller CW vortices at its sides are influenced by the changes in the shear-stress and follow the movement of the central vortex [see Figs. 2 and 5 (Multimedia views)].

For α = 0.45 (supercritical transition from the symmetric to the asymmetric flow), the flow field, prior to the onset of the unsteady flow, contains a monopolar vortex and two co-rotating side vortices (between the co-rotating vortices, there are regions of counter-rotating vorticity). The flow is dominated by three CCW rotating vortices that interact with each other [see Figs. 3 and 6 (Multimedia views)]. The periodic cycles, for α = 0.45, include migration of the central vortex from the center to the side of the yz plane and then a merging event of the central vortex with a co-rotating side vortex that initially formed in the corner of the plane. The initiation of merging between symmetric co-rotating vortices is predicted when the ratio between the vortex radius (R₁) and the separation distance between the vortices (determined between the centers of the cores, b) reaches to a critical value of R₁/b = 0.24–0.3.⁷⁷ However, in the experiments presented here, we do not witness a symmetric vortex merging, since the co-rotating vortices have different intensities and radii. To predict the occurrence of asymmetric vortex merging, a mutuality parameter, which involves the ratio between the strain rate and the vorticity of each vortex, can be used to determine the critical values at which vortex merging will occur.^78,79 However, this is not the focus of the current study.

Experimental results show that it is always the top-left side vortex that intensifies and merges with the central vortex, while numerical simulations show that the top-left and bottom-right side vortices intensify and merge with the central vortex in an alternate manner. Interestingly, this fluctuation pattern resembles a similar periodic flow pattern that was recently reported to occur in a T-channel with two inlets and a single outlet of aspect ratio α = 0.5.²³ In that study, it was reported that at Re = 190, there were two steady co-rotating vortices, with a similar intensity along the channel outlet. At Re_p = 190, close to the start of the channel outlet, a periodic merging of the co-rotating vortices commenced, with new co-rotating vortices being formed at the corners of the channel and migrating toward the center at each cycle of the fluctuation. The process appears to be similar to the one that is seen here for α = 0.45.²³ 

In Fig. 7, we plot St as a function of Re (data are taken from the FFT analysis shown in supplementary material, Figs. S4 and S7). For comparison, we also plot the results from a recent study where the transition to the time-periodic flow was reported to occur in a cross-slot channel with a characteristic length scale of w = 10 mm, α = 1, and Re > 300.²² 

Although all of these flows occur in cross-slot intersections, there are clear differences in the frequencies and the response of the flow field to an increase in the Reynolds number, as found in our work and in the findings from Zhang et al.²² As discussed in Sec. III A, the steady-state flow field is fundamentally different between α = 0.45 and α = 2.4, and the steady-state flow field of α = 1 is also unique.

For the aspect ratio α = 1 (subcritical transition to the unstable flow), the flow field before the onset of the unsteady flow also contains a central vortex and two side vortices. The fluctuations were described as events of continuous merging of two co-rotating Dean vortices and splitting of the central vortex. In the analysis of Zhang et al., they found that the period time of the fluctuation increased with an increase in Re due to a slowdown in the stage of the vortex merging. They reported that the stage in which the flow exhibited a merged central vortex remained nearly constant with a value of ≃0.5 s.²² In their study, the Dean vortices remained weak as Re was increased, which gives a possible explanation for the slowdown in the merging process that was observed. We should also note that their visualization technique, using fluorescent dye from one inlet and non-dyed fluid from the other, does not permit capturing the full complexity of the flow. Using such a method only allows resolving the flow patterns at the interface where the two inflows meet, and it is not possible to resolve side vortices that are formed by the fluid from a single inlet. Therefore, they were not able to visualize the Dean vortices at the symmetric state. The St number that was found through their experiments (0.06 < St < 0.17) is comparable to our results.

The important difference between the steady flow field of α = 1 and the flow field observed for α = 0.45 is that the side vortices seen for α = 1 are counter-rotating in relation to the central vortex (the steady flow field contains CW rotating vortex and two CCW rotating side vortices), while for the case of α = 0.45, the side vortices at the steady state are co-rotating in relation to the central vortex (the steady flow field contains three CCW rotating vortices and two smaller CW rotating vortices between them). This difference strongly affects the merging and splitting dynamics between the vortices in the flow field and can explain the different nature of the periodic fluctuations.

In Fig. 7, the solid black lines are found by fitting the data with Eq. (9). This relation between Re and St was previously used to describe the dimensionless frequency of vortex shedding for flows around circular cylinders with different curvatures, aspect ratios, and angled end plates,^26,80,81

where x₁, x₂, and x₃ are free parameters. The parameters are extracted from the fit and are listed in Table II. For α = 0.45 and α = 1, St as a function of Re follows a near-linear trend.

Periodic vortex flows at the wakes of circular cylinders are characterized by the formation of spanwise vortices, while at intersecting geometries, the flow field consists of streamwise vortices. The quadratic relation in Eq. (9) is used to distinguish between modes of vortex shedding in flows around circular cylinders. Even though vortex properties (i.e., spanwise or streamwise with different core structures) are strongly affected by the flow geometry, we show that the same quadratic relation can also be used for characterizing periodic flows at intersections, particularly for α = 2.4 for the range of Reynolds numbers studied.

Since the nature of the fluctuations is dictated by the flow field at the steady state, the relation between Re and St is expected to be different between geometries. Finding the relation between Re and St can give an initial indication of the modes of fluctuations, controlled by Re and α. Due to experimental limitation, we are not able to study a larger range of Re, but even with the parameter range studied, we are able to characterize our flows and to reveal the similarities between unsteady flows around cylinders and in intersections.

Vortex structure and dynamics of periodic flow fluctuations are studied with two microfluidic cross-slot devices with aspect ratios of α = 2.4 and α = 0.45. Our analysis shows that the vortex in the cross-slot device is similar to a Burgers’ vortex (stretched streamwise vortex) with an additional confinement component that was only previously reported in pancake-like vortices (which are confined in both streamwise and spanwise directions). We find that the aspect ratio changes the properties of streamwise vortices at the steady state. Above Re_c, for α = 0.45, the flow field consists of a central vortex and two, less intense, co-rotating side vortices. Under these conditions, the vorticity across the centerline of the vortex core shows a single Gaussian peak that can be described by a simple pancake-like vortex model. For the higher aspect ratio α = 2.4, the flow field contains a single streamwise vortex (above Re_c) surrounded by a counter-rotating vorticity field. At Re ≳ Re_c, the central vortex has a characteristic single peak in the vorticity across the core centerline (similar to the steady-state flow field for α = 0.45). Further increases in Re while the flow is still steady result in a bubble-type vortex breakdown⁵⁴ and a recirculating flow of the fluid in the upstream direction along the axis of the vortex. This results in a low vorticity center surrounded by a ring of higher intensity vorticity; hence, a double peak in the vorticity profile emerges. In order to describe this vortex, the pancake-like vortex model was modified by the addition of an extra term describing a smaller coaxial vortex core of opposite sign. The modified pancake model was able to accurately capture the complexity of the vorticity distribution at the vortex core.

A further increase in Re up to a critical value, Re_p, results in the onset of a periodic flow cycle. The nature of the periodic cycle depends on the characteristics of the vortical flow field in the steady state. For α = 2.4, upon Re_p, the periodic fluctuations will be pronounced at the vortex core as a collapse and reformation of the core’s structure. Throughout the periodic cycle, the position of the vortex will remain at the center of the geometry. For α = 0.45 upon Re_p, the central vortex will migrate away from the center of the plane (where it was initially formed), where it will merge with one of the co-rotating side vortices in the flow field. The merged vortex will migrate back to the center of the plane and the periodic cycle will start again.

The characteristic frequencies of the fluctuations are found by applying FFT analysis to a velocity time series and are compared to another study that investigated similar fluctuations in the cross-slot geometry with α = 1. We find that the relation between St and Re for α = 2.4 is quadratic, as was previously described for similar stagnation point fluctuations around circular cylinders.^26,80,81 We also find that the relation between St and Re is dependent on the geometry studied, which determines the flow field in the steady state.

The findings of this study show that the frequency of unsteady stagnation point flows can be precisely tuned by slight modifications of Re and α. These results can be directly applied for optimizing applications of inertial microfluidics. Our findings, which are presented in a dimensionless form, are also valuable for improving the control over other stagnation point flow instabilities in various length scales, such as transported fluids within pipes and risers as well as vortex-induced vibrations that may arise in the proximity of obstacles such as bridges, wind turbines, and buildings. Additionally, these results may contribute to increase the efficiency of energy harvesting from vortex-induced vibrations in the proximity of stagnation points, as a source of renewable energy.⁸² 

See the supplementary material for further details of the experiments and simulations referred to in the main text. The figures contain frequency plots and PSD of additional experiments and numerical simulations.

The data that support the findings of this study are available within the article and its supplementary material.

N.B., A.Q.S., and S.J.H. gratefully acknowledge the support of the Okinawa Institute of Science and Technology Graduate University, with subsidy funding from the Cabinet Office, Government of Japan, and funding from the Japan Society for the Promotion of Science (Grant Nos. 17K06173, 17J00412, 18K03958, 18H01135, and 20K14656), and the Joint Research Projects supported by the Japan Society for the Promotion of Science, and the Swiss National Science Foundation. We are grateful to Mr. Kazumi Toda-Peters (Micro/Bio/Nanofluidics Unit, OIST) for help in device fabrication using LightFab. R.J.P. and K.Z. would like to thank the EPSRC for the award of a fellowship under Grant No. EP/M025187/1. N.B. is also grateful to Professor Anke Lindner for helpful discussions and kind hosting at the PMMH lab (ESPCI).