Control of a purely elastic symmetry-breaking flow instability in cross-slot geometries

Abstract

The cross-slot stagnation point flow is one of the benchmark problems in non-Newtonian fluid mechanics as it allows large strains to develop and can therefore be used for extensional rheometry measurements or, once instability arises, as a mixing device. In such a flow, beyond a critical value for which the ratio of elastic force to viscous force is high enough, elasticity can break symmetry even in the absence of significant inertial forces (i.e. creeping flow), which is an unwanted phenomenon if the device is to be used as a rheometer but beneficial from a mixing perspective. In this work, a passive control mechanism is introduced to the cross-slot by adding a cylinder at the geometric centre to replace the ‘free’ stagnation point with ‘pinned’ stagnation points at the surface of the cylinder. In the current modified geometry, effects of the blockage ratio (the ratio of the diameter of the cylinder to the width of the channel), the Weissenberg number (the ratio of elastic forces to viscous forces) and extensibility parameters ($\unicode[STIX]{x1D6FC}$ and $L^{2}$) are investigated in two-dimensional numerical simulations using both the simplified Phan-Thien and Tanner and finitely extensible nonlinear elastic models. It is shown that the blockage ratio for fixed solvent-to-total-viscosity ratio has a stabilizing effect on the associated symmetry-breaking instability. The resulting data show that the suggested modification, although significantly changing the flow distribution in the region near the stagnation point, does not change the nature of the symmetry-breaking instability or, for low blockage ratio, the critical condition for onset. Using both numerical and physical experiments coupled with a supporting theoretical analysis, we conclude that this instability cannot therefore be solely related to the extensional flow near the stagnation point but it is more likely related to streamline curvature and the high deformation rates towards the corners, i.e. a classic ‘curved streamlines’ purely elastic instability. Our work also suggests that the proposed geometric modification can be an effective approach for enabling higher flow rates to be achieved whilst retaining steady symmetric flow.

1 Introduction

Extensional flow occurs when a fluid experiences a deformation in the streamwise direction and can be observed in many different situations such as flow passing through contractions (Afonso et al. 2011) and expansions (Pinho, Oliveira & Miranda 2003; Alves & Poole 2007) and at stagnation point flows like flow through intersections (such as T- shaped (Soulages et al. 2009) or cross-slot (Arratia et al. 2006) geometries) or flow passing an obstacle (Bisgaard 1983; Walters & Tanner 1992; Alves, Pinho & Oliveira 2001; McKinley 2002). Many industrial processes deal with highly extensional flows of polymeric fluids. In most of these cases, a purely elastic instability occurs that is absent for its equivalent Newtonian creeping flow (Poole et al. 2007a ; Zilz et al. 2014). During the last few decades, significant attention has been placed on such flows to characterize the physics behind these types of instabilities especially in shear flows (Larson, Shaqfeh & Muller 1990; McKinley, Pakdel & Öztekin 1996); however, understanding of this type of elastic instability in extensional-dominated flows is not as advanced as for shear-dominated flows (Haward, McKinley & Shen 2016).

The cross-slot geometry is a common flow geometry that has been utilized for generating controllable planar elongational flow fields and for studying the stretching dynamics of polymers (Haward et al. 2012b , 2016). This geometry consists of four bisecting rectangular channels with two sets of opposite inlets and outlets. These opposing inlets and outlets produce a flow field with a free stagnation point located at the centre of the geometry. At this point the velocity field is zero and a finite gradient of velocity in the streamwise direction appears. In principle, due to the zero velocity field at this point, a fluid element is trapped for an infinite time generating a significant strain and potentially enabling ‘steady-state’ extentional flow kinematics to be realized. To obtain an ideal planar elongational flow field, Haward et al. (2012b ) suggested optimizing the standard shape of the cross-slot geometry using a numerical approach and termed the resulting shape the optimized shape cross-slot elongation rheometry (OSCER) geometry. A similar type of stagnation point flow can be observed when fluid is passing an object such as a cylinder or a sphere (Walters & Tanner 1992). If the obstacle is a solid object the stagnation point is pinned at the surface of the geometry and cannot move, while for example in the case of falling/rising drops, the stagnation point is located at the surface of the moving drop and in principle is free to move and change the shape of the droplet (Bisgaard 1983; McKinley 2002; Davoodi & Norouzi 2016). As a consequence of the no-slip condition and mass continuity, the local velocity gradient and the velocity field are zero for a pinned stagnation point while for a free stagnation point the strain rate can have finite non-zero values. In both types of flows, beyond a critical value of strain rate a symmetry breaking of the flow distribution can be observed for both the pinned (Hulsen, Fattal & Kupferman 2005) and free (Bisgaard 1983) stagnation point flows. Soulages et al. (2009) studied two different geometries to investigate the kinematic differences between a pinned stagnation point flow at a wall and a free stagnation point flow by adding a recirculating cavity opposite to the outlet arm of a T-shaped channel. They showed that the critical value of flow rate for the onset of an unsteady three-dimensional instability is delayed to higher values of the Weissenberg number for the free stagnation point flow geometry with a cavity compared to the pinned case without a cavity. Also, a new type of steady asymmetric instability was reported in this modified geometry that was suppressed for the cases in which the stagnation point was pinned at the wall.

Early research conducted by Gardner et al. (1982) was the first to report that a steady flow asymmetry can occur for viscoelastic flows in cross-slot geometries. In this geometry, although nominally extensional dominated, fluid particles passing through the cross-section between the corner and the stagnation point experience both significant shear flow near these re-entrant corners (Dean & Montagnon 1949; Moffatt 1964; Davies & Devlin 1993; Hinch 1993) and elongational-dominated flow near the stagnation point (Öztekin, Alakus & McKinley 1997; Haward et al. 2016). The combination of this complex deformation with the nonlinear elastic stresses for viscoelastic materials can enable disturbances to grow and trigger a ‘purely elastic instability’. Although firstly observed by Gardner et al. (1982), it was not until Arratia et al. (2006) that this effect was clearly associated with a purely elastic instability and suitably characterized. Supporting numerical simulations for this phenomenon were presented by Poole, Alves & Oliveira (2007b ), where it was shown that once the instability appears, the shape of the velocity profile along a line between the corner of the cross-slot geometry and the free stagnation point changes from convex into concave. These types of instability are generally triggered when a combination of the normal stress in the streamwise direction is coupled with streamline curvature, although in the cross-slot geometry this is debated (Wilson 2012) and something we will directly address in the current article. A well-known dimensionless parameter which rationalizes these ‘curved streamlines’ instabilities is the $M$ parameter introduced by McKinley et al. (1996) (often referred to as the ‘Pakdel–McKinley’ criteria). This parameter can be considered as the viscoelastic complement of the Görtler number (Saric 1994) and is defined as

(1.1) $$\begin{eqnarray}\displaystyle M=\sqrt{{\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}\tilde{U}}{\tilde{\Re }}}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D70F}}_{ss}}{\tilde{\unicode[STIX]{x1D702}}_{0}\tilde{\dot{\unicode[STIX]{x1D6FE}}}}}}, & & \displaystyle\end{eqnarray}$$

where $\tilde{U}$ is the magnitude of the local velocity, $\tilde{\Re }$ is the local radius of curvature of a streamline, $\tilde{\unicode[STIX]{x1D70F}}_{ss}$ is the normal stress in the streamwise direction, $\tilde{\unicode[STIX]{x1D702}}_{0}$ is the zero-shear-rate viscosity of the fluid and $\tilde{\dot{\unicode[STIX]{x1D6FE}}}$ is the magnitude of the shear rate. Throughout this paper we will indicate dimensional variables using a tilde. In (1.1), the first term on the right-hand side can be referred to as a local Deborah number, showing the ratio of the relaxation time of the fluid to the time a disturbance takes to travel along a streamline. As this ratio increases, the chance that a disturbance may grow and lead to instability increases. The second term on the right-hand side of (1.1) is added to scale properly the effect of the normal stress in the streamwise direction with a reference stress scale. This term is generally of the same order of magnitude as a local Weissenberg number. Recently, Cruz et al. (2016) have plotted the spatial variation of this parameter in the cross-slot geometry for both upper-convected Maxwell and simplified Phan-Thien–Tanner (sPTT) fluids. Previous studies for flow around sharp re-entrant corners for both Newtonian (Dean & Montagnon 1949; Moffatt 1964) and upper-convected Maxwell (Davies & Devlin 1993; Hinch 1993) fluids suggest that although the velocity field at the corner is zero, the velocity gradient, and consequently the magnitude of the stress tensor, are singular. At such corners the combination of high shear rate along with high curvature of the streamlines can thus provide a suitable mechanism for a disturbance to grow and trigger an instability. The numerical simulations of Cruz et al. (2016) for a geometry with sharp square re-entrant corners (i.e. a standard-shape cross-slot) suggest that instability should be triggered at a region near the corner of the geometry. These results support the earlier work by Rocha et al. (2009) who showed the instability is delayed to a higher value of the Weissenberg number once the sharp square corner is replaced with a strongly rounded corner. In a following work, Haward et al. (2016) using the OSCER geometry showed that once the ideal planar elongation flow field is obtained, although the velocity exhibits a small value at the stagnation point, the curvature of the streamline exhibits a large value near this region and the maximum value of the $M$ parameter appears in the vicinity of the stagnation point. Despite these results, it still remains an open question in the literature as to whether instabilities in the cross-slot are driven via curvature and high deformation rate near the stagnation point or from a region closer to the corner, or from simply the extensional flow at the stagnation point itself (Wilson 2012; Kalb, Villasmil & Cromer 2017).

In the current work, using a series of numerical simulations supported by experimental and analytical tools, we study and modify the standard-shape cross-slot geometry by replacing the free stagnation point flow with pinned stagnation points by adding a cylinder at the geometric centre of the domain. Using this modification, the curvature of streamlines close to the geometric centre will be affected by the presence of the cylinder and will exhibit an opposite sign in comparison to the streamlines near the corners of the geometry. A change in the relative size of the cylinder (i.e. changing the ‘blockage ratio’ parameter) will allow us to control the local value of streamline curvature along with the average value of the velocity passing through the gap between the corner of the cross-slot and the cylinder. We will show that the introduced modification acts as a passive control mechanism that can be used to delay the critical flow rate in which the instability is triggered. Also, due to conservation of mass and the no-slip boundary condition, once a cylinder is added, the finite non-zero value of the strain rate at the stagnation point seen for the standard-shape cross-slot geometry is replaced by stagnation points where the strain rates are identically zero. Considering the above, one can state that the addition of the cylinder at the geometrical centre of the cross-slot geometry significantly changes the flow distribution near the stagnation point while the flow distribution close to the corners will stay relatively unchanged (at least for small cylinders).

The structure of the paper is as follows. In § 2, the governing equations and the employed numerical procedure are reviewed. In § 3, the experimental set-up and the employed protocol are explained. Following this, in § 4, a discussion of the results obtained in this study is presented followed by conclusions.

2 Numerical procedure

The governing equations for this problem are conservation of mass, assuming incompressibility, and momentum:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D70C}}\left({\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}}{\unicode[STIX]{x2202}\tilde{t}}}+\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}\right)=-\unicode[STIX]{x1D735}\tilde{p}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D749}}+\tilde{\unicode[STIX]{x1D702}}_{s}\unicode[STIX]{x1D6FB}^{2}\tilde{\boldsymbol{u}}, & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D70C}}$ is the density of the fluid, $\tilde{\boldsymbol{u}}$ is the velocity vector, $\tilde{\unicode[STIX]{x1D702}}_{s}$ is the solvent viscosity and $\tilde{\unicode[STIX]{x1D749}}$ is the extra-stress tensor containing the polymeric contribution to the stress. In this work two different constitutive equations, namely the sPTT model (Phan-Thien & Tanner 1977) and the finitely extensible nonlinear elastic (FENE-P) model (Bird, Dotson & Johnson 1980), have been used to study the effect of different constitutive equations.

The extra-stress tensor using the sPTT model may be calculated as follows:

(2.3) $$\begin{eqnarray}\displaystyle f_{1}\tilde{\unicode[STIX]{x1D749}}+\tilde{\unicode[STIX]{x1D706}}\stackrel{\unicode[STIX]{x1D735}}{\tilde{\unicode[STIX]{x1D749}}}=\tilde{\unicode[STIX]{x1D702}}_{p}(\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}+\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}^{\text{T}}), & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D706}}$ is the relaxation time of the fluid and $\tilde{\unicode[STIX]{x1D702}}_{p}$ is the polymeric contribution to the viscosity. The upper-convective derivative of the extra-stress tensor, $\stackrel{\unicode[STIX]{x1D735}}{\tilde{\unicode[STIX]{x1D749}}}$ , is defined as

(2.4) $$\begin{eqnarray}\displaystyle \stackrel{\unicode[STIX]{x1D735}}{\tilde{\unicode[STIX]{x1D749}}}={\displaystyle \frac{\text{D}}{\text{D}\tilde{t}}}(\tilde{\unicode[STIX]{x1D749}})-(\tilde{\unicode[STIX]{x1D749}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}+\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}^{\text{T}}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D749}}), & & \displaystyle\end{eqnarray}$$

where the material derivative of an arbitrary matrix $\tilde{\boldsymbol{A}}$ is defined as $(\text{D}/\text{D}\tilde{t})\tilde{\boldsymbol{A}}=(\unicode[STIX]{x2202}\tilde{\boldsymbol{A}}/\unicode[STIX]{x2202}\tilde{t}+\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{\boldsymbol{A}})$ . The $f_{1}$ function for the linear sPTT model is defined as follows:

(2.5) $$\begin{eqnarray}\displaystyle f_{1}=1+\unicode[STIX]{x1D6FC}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}}{\tilde{\unicode[STIX]{x1D702}_{p}}}}Tr(\tilde{\unicode[STIX]{x1D749}}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is the the extensibility parameter, and, in the limiting case that $\unicode[STIX]{x1D6FC}=0$ , the sPTT constitutive equation reduces to the Oldroyd-B model. This constitutive equation is derived from network theory by Phan-Thien & Tanner (1977) and is a suitable model for simulation of shear-thinning polymeric fluids (Bird, Armstrong & Hassager 1987).

One other constitutive equation which has frequently been used to simulate the shear-thinning behaviour of polymeric materials is the FENE-P model (Bird et al. 1980) which may be presented in the following form:

(2.6) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D749}}+\tilde{\unicode[STIX]{x1D706}}\stackrel{\unicode[STIX]{x1D735}}{\left({\displaystyle \frac{\tilde{\unicode[STIX]{x1D749}}}{f_{2}}}\right)}={\displaystyle \frac{a\tilde{\unicode[STIX]{x1D702}}_{p}}{f_{2}}}(\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}+\unicode[STIX]{x1D735}\tilde{\boldsymbol{u}}^{\text{T}})-a\tilde{\unicode[STIX]{x1D702}}_{p}\left({\displaystyle \frac{\text{D}}{\text{D}\tilde{t}}}\left({\displaystyle \frac{1}{f_{2}}}\right)\right)\boldsymbol{I}, & & \displaystyle\end{eqnarray}$$

where $a=L^{2}/(L^{2}-3)$ , $\boldsymbol{I}$ is the identity tensor and $L^{2}$ is also called the extensibility parameter. The $f_{2}$ function for the FENE-P model is defined as

(2.7) $$\begin{eqnarray}\displaystyle f_{2}={\displaystyle \frac{L^{2}+{\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}}{a\tilde{\unicode[STIX]{x1D702}}_{p}}}Tr(\tilde{\unicode[STIX]{x1D749}})}{L^{2}-3}}. & & \displaystyle\end{eqnarray}$$

Here we use the rheoFoam solver in the OpenFOAM platform which was previously introduced by Pimenta & Alves (2017). In this solver, instead of dealing with large values of stress components in the global $\tilde{x}$ , ${\tilde{y}}$ and $\tilde{z}$ coordinate system, the logarithms of the eigenvalues of the stress tensor in a local coordinate system consisting of the eigenvectors of the stress tensor (i.e. principal axis) are calculated and solved. More detail about the employed approach can be obtained from Afonso, Pinho & Alves (2012) and Fattal & Kupferman (2004).

2.1 Problem definition

A schematic of the geometry and the employed coordinate systems are shown in figure 1. The length of the inlet/outlet arms is set to be 15 times the width of the channel in order to allow fully developed conditions at the cross-slot, which we confirmed in our numerical simulation to be sufficient. In the numerical procedure, a constant bulk velocity $\tilde{U}_{B}$ at the inlets and Neumann boundary condition at the outlets were applied. At the walls, no-slip boundary conditions were imposed and the values of the extra stress components were linearly extrapolated using the method introduced in Pimenta & Alves (2017). In order to better understand effects of the proposed geometry modification on the cross-slot symmetry-breaking instability, besides the no-slip boundary condition, a number of additional simulations using a complete slip boundary condition at the cylinder wall were also carried out (i.e. normal velocity component is fixed to zero while the tangential component is set to have zero gradient). No finite disturbances are introduced in the simulations to induce the onset of symmetry-breaking instability. Instead, the instability is naturally triggered from accumulation of numerical error via machine-level precision.

Figure 1. Schematic of the cross-slot geometry with a cylinder and the coordinate systems used in the problem. Not to scale. The angle $\unicode[STIX]{x1D703}$ is measured with respect to the diagonal line $|x|=|y|$ .

2.2 Non-dimensionalization

It is convenient to use dimensionless parameters in this problem. The relationships between dimensional and dimensionless parameters are

(2.8) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}x={\displaystyle \frac{\tilde{x}}{\tilde{W}}},\quad y={\displaystyle \frac{{\tilde{y}}}{\tilde{W}}},\quad r={\displaystyle \frac{\tilde{r}}{\tilde{W}}},\quad \unicode[STIX]{x1D6F7}={\displaystyle \frac{\tilde{D}}{\tilde{W}}},\quad \boldsymbol{U}={\displaystyle \frac{\tilde{\boldsymbol{U}}}{\tilde{U}_{B}}},\quad Re={\displaystyle \frac{\tilde{\unicode[STIX]{x1D70C}}\tilde{U}_{B}\tilde{W}}{\tilde{\unicode[STIX]{x1D702}}_{t}}},\\[12.0pt] Wi={\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}\tilde{U}_{B}}{\tilde{W}}},\quad \unicode[STIX]{x1D6FD}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{s}}{\tilde{\unicode[STIX]{x1D702}}_{t}}},\quad N_{1}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D70F}}_{xx}-\tilde{\unicode[STIX]{x1D70F}}_{yy}}{\tilde{\unicode[STIX]{x1D702}}_{t}\tilde{U}_{B}/\tilde{W}}},\quad \dot{\unicode[STIX]{x1D716}}={\displaystyle \frac{\tilde{\dot{\unicode[STIX]{x1D716}}}}{\tilde{U}_{B}/\tilde{W}}},\quad \unicode[STIX]{x1D734}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D734}}}{\tilde{U}_{B}/\tilde{W}}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

were $\tilde{x}$ , ${\tilde{y}}$ and $\tilde{r}$ are the variables related to the corresponding rectangular and polar coordinate systems, $\tilde{W}$ is the width of the channel, $\tilde{D}$ is the diameter of the cylinder, $\unicode[STIX]{x1D6F7}$ is the blockage ratio parameter which may change between zero (the standard cross-slot geometry with no cylinder) and $\sqrt{2}$ (fully blocked cross-slot geometry), $\tilde{\boldsymbol{U}}$ is the velocity vector, $\tilde{U}_{B}$ is the imposed bulk velocity at the inlet arms, $Re$ is the Reynolds number which is set to be 0.01 for all simulations in order to model creeping flow, $Wi$ is the Weissenberg number, $\unicode[STIX]{x1D6FD}$ is the solvent-to-total-viscosity ratio, $\tilde{\unicode[STIX]{x1D702}_{t}}$ is the total viscosity (i.e. $\tilde{\unicode[STIX]{x1D702}}_{t}=\tilde{\unicode[STIX]{x1D702}}_{s}+\tilde{\unicode[STIX]{x1D702}}_{p}$ ), $N_{1}$ is the non-dimensional first normal-stress difference, $\tilde{\dot{\unicode[STIX]{x1D716}}}$ is the strain rate (either $\tilde{\dot{\unicode[STIX]{x1D716}}}_{xx}$ or $\tilde{\dot{\unicode[STIX]{x1D716}}}_{yy}$ ) and $\tilde{\unicode[STIX]{x1D734}}$ is the vorticity tensor.

2.3 Mesh dependency study

In this section, a number of representative results analysing the effect of mesh on the two-dimensional (2-D) flow distribution are presented to give an overview of numerical accuracy. The block-structured mesh generation in OpenFOAM has required us to divide the flow domain into 16 smaller sub-domain blocks. In figure 2, a schematic definition of different blocks used in our mesh-generation steps for a nominal blockage ratio of $\unicode[STIX]{x1D6F7}=0.55$ is illustrated. Here, two different meshes were used and the characteristics of these meshes are presented in table 1. The cell growth rate is defined based on the ratio of the first cell size to the last cell size in a specific direction. Note that, owing to the symmetry of the domain, only mesh information for a quarter of the geometry is presented but the full geometry is used in all simulations.

Figure 2. Schematic representation of the geometry and definition of the blocks used to generate the mesh.

Table 1. Main characteristics of meshes in different blocks (NR $\times$ NT ^a (GRR, GRT) ^b problems).

a NR, NT, number of cells in radial and tangential directions; NCV, total number of control volumes.

b GRR, GRT, cell growth rate in radial and tangential directions.

Figure 3. Distribution of velocity magnitude at the entrance region of the outlet arms (i.e. along the line $y=\pm 0.5$ ) for different meshes for Newtonian creeping flow. Results for M1 and M2 are shown in blue and red, respectively.

In figure 3, the effect of mesh refinement on the flow distribution of creeping Newtonian flow is presented (i.e. $Re=0.01$ ). In this figure, the distribution of velocity magnitude at the entrance region of the outlet arms is presented for the meshes introduced in table 1. The results show that the mean average errors between M1 and M2 are less than $1.3\,\%$ . Therefore, the remaining results presented here were obtained using the M1 mesh, except for simulations corresponding to $\unicode[STIX]{x1D6F7}=0.10$ in which meshes similar in density to M2 were used. For the standard cross-slot geometry a mesh similar to that of Cruz et al. (2016) consisting of five blocks was used. The total number of cells for the standard cross-slot geometry is 13 005.

3 Experiments

3.1 Experimental set-up

A schematic diagram of the experimental rig is shown in figure 4. Flow through the microfluidic system was driven using two identical individually controlled high-precision syringe pumps (PHD Ultra Harvard Apparatus). One of the pumps drives fluid into the two opposed inlets, while the other one withdraws fluid simultaneously from the two outlets of the device (all at equal volumetric flow rates). According to the manufacturer, the mass flow rate ( $\tilde{{\dot{m}}}$ ) certified accuracy is $0.35\,\%$ at the lowest pulsation-free delivering rate, which set our lowest flow rate.

Figure 4. (a) A schematic illustrating the experimental apparatus of a microfluidic cross-slot device allowing for direct observation of the $x{-}y$ plane. The origin is placed at the geometric centre of the device. The rig is mounted on an inverted microscope fitted with a filter cube. A pulsed Nd:YAG laser is used to excite the dyed fluid, and a CCD camera enables the instability formation to be captured (the illustration is not to scale; the two-camera field of view (FOV) to be captured is $4~\text{mm}\times 4~\text{mm}$ ). (b) The prototype microfluidic device rigs for (I) standard cross-slot geometry and (II) cross-slot with cylinder geometry. The channels were micro-machined in brass and encased in polyoxymethylene. (c) Photograph illustrating the assembled experimental rig set-up.

For the purpose of this study, two geometries, namely one standard cross-slot device and one modified cross-slot geometry with $\unicode[STIX]{x1D6F7}=0.55$ , were designed. These geometries were micro-machined into two pieces of brass using CNC machining with a rectangular cross-section channel. To best approximate a 2-D channel, the channel height is selected to be twice its width, $\tilde{H}=2000\pm 10~(\unicode[STIX]{x03BC}\text{m})$ and $\tilde{W}=1000\pm 10~(\unicode[STIX]{x03BC}\text{m})$ , respectively. The device was encased in polyoxymethylene, an insulating material also known as ACETAL. The cross-sectional dimensions of the cross-slot flow device were quantified using a Nikon EPIPHOT TME inverted microscope with $\times$  100 magnification, $470~\text{pixels}=1000~(\unicode[STIX]{x03BC}\text{m})$ . The combined length of the channel inlet and outlet is $80\tilde{W}=80\pm 0.1~(\text{mm})$ to ensure fully developed flow at the central region of the geometry for all flow rates studied. The cross-slot geometry was enclosed by a 6.5 (mm) thick upper wall fabricated from borosilicate glass to maintain sealing while still allowing the flow structure to be visualized.

Rhodamine-B (Acros Organics) was chosen as a fluorescent dye to capture the flow patterns in the cross-slot micro-geometry (Ross, Gaitan & Locascio 2001; Huang et al. 2013, 2014). The dyed fluid was prepared by dissolving 30 ppm of Rhodamine-B in half of the working fluid. The optical set-up was composed of an inverted microscope (DMI, Leica Microsystems GmbH) fitted with an appropriate filter cube (excitation BP 530–545 (nm), dichroic mirror 565 (nm), barrier filter 610–675 (nm), Leica Microsystems GmbH), a pulsed Nd:YAG laser (Solo-PIV III laser, wavelength 532 (nm), New Wave Research) used to excite the fluorescent Rhodamine-B dye with illumination and a CCD (charge couple device) camera.

3.2 Working fluids

The Newtonian working fluid used in the experiments is a mixture of glycerine (relative density 1.26, ReAgent Chemical Services) and distilled water with a nominal concentration of 70 % glycerine (by weight). The density ( $\tilde{\unicode[STIX]{x1D70C}}$ ) and viscosity ( $\tilde{\unicode[STIX]{x1D702}}_{s}$ ) of the mixed solution were measured at $\tilde{T}=20~(\text{}^{\circ }\text{C})$ . A density meter (Anton Paar DMA 35N) with a quoted precision of $0.001~(\text{g}~\text{cm}^{-3})$ was used for quantifying the fluid density. A controlled-stress torsional rheometer (Anton Paar MCR302) was utilized to measure the fluid viscosity using a 60 (mm) and $1^{\circ }$ cone with a shear rate range from $1\times 10^{-3}(\text{s}^{-1})$ to $200~(\text{s}^{-1})$ . The measured values of the density ( $\tilde{\unicode[STIX]{x1D70C}}$ ) and viscosity ( $\tilde{\unicode[STIX]{x1D702}}_{s}$ ) are $1184~(\text{kg}~\text{m}^{-3})$ and $35.50~(\text{mPa}~\text{s})$ , respectively.

Figure 5. Analytical fit based on the sPTT model ( $\unicode[STIX]{x1D6FC}=0.63,\tilde{\unicode[STIX]{x1D702}_{0}}=3360~(\text{mPa}~\text{s}),\tilde{\unicode[STIX]{x1D702}}_{\infty }=35.5~(\text{mPa}~\text{s})$ ) and the experimental shear viscosity data of 190 ppm PAA in 70 : 30 glycerine–water solution versus the shear rate at $\tilde{T}=20~(\text{}^{\circ }\text{C})$ .

To elucidate the role of elasticity, a mixture of 70 % glycerine and distilled water (by weight) and 190 ppm polyacrylamide (PAA; $M_{\text{w}}=18\times 10^{6}~(\text{g}~\text{mol}^{-1})$ , Sigma-Aldrich) was used. Initially, the required amount of PAA is added to distilled water, and the resultant liquid is slowly mixed using a magnetic stirrer. In another container, the necessary amount of glycerine is also added to distilled water and mixed using a magnetic stirrer. After 24 h, the result is two clear, transparent and colourless solutions. The PAA and glycerine solutions are then combined in a single container, and the resultant polymer solution is thus left to mix gently in a magnetic stirrer (24 h). All rheology measurements were conducted at $\tilde{T}=20\pm 0.1~(\text{}^{\circ }\text{C})$ . The working fluid characterization includes measuring viscosity, density and relaxation time. Figure 5 displays the experimental data for viscosity of the viscoelastic fluid as a function of applied shear rate and an analytical fit based on a single-mode sPTT model with $\unicode[STIX]{x1D6FC}=0.63,\tilde{\unicode[STIX]{x1D702}}_{0}=3.36~(\text{Pa}~\text{s})$ and $\tilde{\unicode[STIX]{x1D702}}_{\infty }=35.5~(\text{mPa}~\text{s})$ . The relaxation time of the fluid was measured using small-amplitude oscillatory shear. Analysing the storage modulus ( $\tilde{G}^{\prime }$ ) and the loss modulus ( $\tilde{G}^{\prime \prime }$ ) using a Maxwell model is a standard approach for estimating the relaxation time of viscoelastic fluids (Bird et al. 1987). Here, a four-mode generalized Maxwell model is used to estimate the average relaxation time. The response of the generalized Maxwell model for the oscillatory test is as follows:

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{G}^{\prime }=\mathop{\sum }_{i=1}^{4}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{i}\tilde{\unicode[STIX]{x1D706}}_{i}\tilde{\unicode[STIX]{x1D714}}^{2}}{1+\tilde{\unicode[STIX]{x1D706}}_{i}^{2}\tilde{\unicode[STIX]{x1D714}}^{2}}}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{G}^{\prime \prime }=\mathop{\sum }_{i=1}^{4}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{i}\tilde{\unicode[STIX]{x1D714}}}{1+\tilde{\unicode[STIX]{x1D706}}_{i}^{2}\tilde{\unicode[STIX]{x1D714}}^{2}}}, & \displaystyle\end{eqnarray}$$

where the unknown $\tilde{\unicode[STIX]{x1D706}}_{i}$ and $\tilde{\unicode[STIX]{x1D702}}_{i}$ values of the model are determined by minimizing the deviation between the model and experimental data of the frequency sweep test subject to the constraints that the total viscosity is equal to $3.36~(\text{Pa}~\text{s})$ and that of the solvent mode is $35.5~(\text{mPa}~\text{s})$ . Finally, the average relaxation time of our working fluid may be estimated as

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\tilde{\unicode[STIX]{x1D706}}}={\displaystyle \frac{\displaystyle \mathop{\sum }_{i=1}^{4}\tilde{\unicode[STIX]{x1D702}}_{i}\tilde{\unicode[STIX]{x1D706}}_{i}}{\displaystyle \mathop{\sum }_{i=1}^{4}\tilde{\unicode[STIX]{x1D702}}_{i}}}. & \displaystyle\end{eqnarray}$$

As for the shear viscosity, the Anton Paar MCR302 is used to measure the small-amplitude oscillatory shear properties of the working fluid. Both the amplitude and frequency sweep tests are carried out at $T=20~(\text{}^{\circ }\text{C})$ , using a 60 (mm) and $1^{\circ }$ cone. Figure 6 shows the results of the frequency sweep test in the range of 0.01 to $4~(\text{rad}~\text{s}^{-1})$ at $10\,\%$ strain amplitude (which we confirmed was in the linear limit). Our experimental results suggest an average relaxation time of $\bar{\tilde{\unicode[STIX]{x1D706}}}=45.05~(\text{s})$ .

Figure 6. Variation of the storage and loss modulus versus the frequency for the solution of 190 ppm PAA in 70 : 30 glycerine–water at $\tilde{T}=20~(\text{}^{\circ }\text{C})$ .

3.3 Experimental protocol

Experiments are conducted over a range of $0.1<\bar{Wi}<200$ by programming the syringe pumps to perform ramps in $\bar{Wi}$ with small step increases in volumetric flow rate ( $0.001<\dot{Q}<0.2~(\text{nl}~\text{s}^{-1})$ ). Here, $\bar{Wi}$ is defined based on the average relaxation time of the fluid obtained from (3.3). The camera is synchronized with the laser at a repetition rate of 8.875 (Hz). All images shown in this work were captured from the top view using a $\times$  8 microscope objective (Leica Microsystems GmbH) and a camera field of view of $4~\text{mm}\times 4~\text{mm}$ which thoroughly covers the central region of the device containing the field of interest. The flow visualization technique applied consists of pumping dyed fluid (Rhodamine-B) from one of the inlets while undyed fluid is pumped in the other inlet; both inlets are kept at the same flow rate using the same syringe pump and two identical syringes. Fluid motion starts from rest and at least 50 images were captured for each steady state. Due to the flow visualization technique used the images represent an integral measurement across a significant portion of the channel depth. To ensure the flow distribution has reached steady state, after changing the flow rate we wait approximately 30–40 minutes at each step, i.e. more than 45 fluid relaxation times.

4 Results and discussion

In this section, results related to the influence of the geometric modification on the flow field and on the associated steady symmetry-breaking instability in the cross-slot geometry are presented. We use a series of 2-D numerical simulations, experimental results and a supporting approximate analytical solution to determine the $M$ parameter (1.1) to show that this instability cannot be solely related to the extensional flow near the stagnation point and is more likely related to streamline curvature and the high deformation rates nearer the re-entrant corners.

Figure 7. Distribution of the flow-type parameter and the non-dimensional $xy$ component of the vorticity tensor ( $\unicode[STIX]{x1D6FA}_{xy}$ ) for Newtonian and viscoelastic fluids with $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ , in a standard-shape cross-slot geometry ( $\unicode[STIX]{x1D6F7}=0$ ). Streamlines are superimposed in all cases.

Figure 8. Comparison of non-dimensional velocity distribution along (a) the vertical line at $x=-0.5$ , (b) the horizontal line at $y=0.5$ and (c) the horizontal line at $y=1$ between Newtonian and viscoelastic fluids with $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ , in a standard-shaped cross-slot geometry ( $\unicode[STIX]{x1D6F7}=0$ ).

In the standard cross-slot geometry, fluid particles experience a complex deformation due to the existence of shear-dominated and elongational-dominated flows near the corners and stagnation point, respectively. A suitable parameter that can be used to visualize different types of flow is the flow-type parameter $\unicode[STIX]{x1D709}$ (Lee et al. 2007):

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}={\displaystyle \frac{||\tilde{\boldsymbol{D}}||-||\tilde{\unicode[STIX]{x1D734}}||}{||\tilde{\boldsymbol{D}}||+||\tilde{\unicode[STIX]{x1D734}}||}}, & \displaystyle\end{eqnarray}$$

where $||\tilde{\boldsymbol{D}}||=((\tilde{\boldsymbol{D}}:\tilde{\boldsymbol{D}})/2)^{1/2}$ and $||\tilde{\unicode[STIX]{x1D734}}||=((\tilde{\unicode[STIX]{x1D734}}:\tilde{\unicode[STIX]{x1D734}})/2)^{1/2}$ are magnitudes of the rate of deformation (i.e. $\tilde{\boldsymbol{D}}=(\unicode[STIX]{x1D735}\tilde{\boldsymbol{U}}+\unicode[STIX]{x1D735}\tilde{\boldsymbol{U}}^{\text{T}})/2$ ) and vorticity tensors (i.e. $\tilde{\unicode[STIX]{x1D734}}=(\unicode[STIX]{x1D735}\tilde{\boldsymbol{U}}-\unicode[STIX]{x1D735}\tilde{\boldsymbol{U}}^{\text{T}})/2$ ), respectively. The $\unicode[STIX]{x1D709}$ parameter may vary within the range $[-1,1]$ in which $\unicode[STIX]{x1D709}=-1$ characterizes a solid-like rotational flow, $\unicode[STIX]{x1D709}=1$ a pure extensional flow and $\unicode[STIX]{x1D709}=0$ a simple shear flow. Results presented in the left-hand column of figure 7 show the effect of elastic stress on the streamline distribution superimposed on the flow-type parameter for small values of $Wi$ before the onset of the symmetry-breaking purely elastic instability. In all figures, the inlets are located on the left- and right-hand sides while the outlets are at the top and bottom. Due to the existence of these opposite inlets and outlets, a point with zero velocity appears at the centre of the geometry (a stagnation point) resulting in a planar elongational flow field in this area. As shown in figure 7, in the creeping flow of Newtonian fluids, this extensional-dominated flow appears in the shape of four strands stretched along the centrelines of the inlet/outlet arms. These strands are located at the mid-distance between the two walls where, due to symmetry, the shear rate is zero (see figure 8), and due to a non-zero value of the streamwise gradient of velocity, a purely elongational flow is observed. This effect is also highlighted by plotting the $xy$ component of the vorticity tensor ( $\unicode[STIX]{x1D6FA}_{xy}$ ) as shown in the right-hand column of figure 7.

On increasing the Weissenberg number up to a certain limit (here, $Wi=0.17$ for $\unicode[STIX]{x1D6FD}=1/9$ and $\unicode[STIX]{x1D6FC}=0.02$ ), the velocity distribution in the fully developed inlet arms exhibits a flatter distribution in comparison to the Newtonian flow field, due to shear-thinning, and consequently the elongational-dominated flow expands to a wider region. As shown in figure 8, on increasing the value of the Weissenberg number beyond this limit, three shear-free locations appear in the velocity profile, leading to three strands of elongational-dominated fields at the outlet arms (along the line $y=\pm 0.5$ ). As the flow redevelops we expect the flow further downstream to eventually exhibit its approximately parabolic fully developed distribution (i.e. with only one shear-free point), so these strands naturally join together as we move further downstream through the outlet arms. The results in figures 7 and 8 show that the effect of elasticity on the flow field, even prior to any purely elastic instability, is significant and changes the flow type downstream of the stagnation point in a quite fundamental way.

Figure 9. Distribution of the flow-type parameter for Newtonian and viscoelastic fluids with $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ .

Previous experimental studies (Gardner et al. 1982; Arratia et al. 2006) conducted with the standard cross-slot geometry have shown that beyond a critical value of $Wi$ , elasticity can break symmetry even in the creeping-flow regime. This change of flow distribution can be related to the existence of multiple solutions to the hydrodynamic problem due to the nonlinear nature of the elastic force in viscoelastic fluids. In other words, due to the existence of the nonlinear upper-convective derivative operator (2.4) in the governing equations, the problem may exhibit more than one solution leading to bifurcation of solutions after a critical value of Weissenberg number. The modifying effects of adding a cylinder to the flow field and the associated symmetry-breaking instability are presented in figure 9 (for a nominal case of $\unicode[STIX]{x1D6F7}=0.50$ ). The results suggest that once the cylinder is added, although the flow distribution in the geometric central region changes significantly, i.e. the shear-free elongational-dominated regime with a single stagnation point is replaced with a shear flow around the cylinder wall containing four stagnation points stretched towards the inlet/outlet arms, the steady symmetry-breaking nature of the instability essentially does not change.

Figure 10. Effect of geometry modification on streamlines superimposed onto contour plots of non-dimensional first normal-stress difference for $Wi=0.7,\unicode[STIX]{x1D6FC}=0.02$ , $\unicode[STIX]{x1D6FD}=1/9$ . At $\unicode[STIX]{x1D6F7}=0.60$ and beyond, symmetry is returned to the flow.

By increasing the blockage ratio, the average velocity passing through the central region of the cross-slot geometry (i.e. through $\tilde{L}$ in figure 1) increases which consequently increases the local shear rate and hence the first normal-stress difference in this region (see figure 10). On the other hand, by changing the size of the cylinder, the local values of streamline curvature are also changed. These two issues suggest that the proposed geometry modification may be used as an effective approach for controlling the kinematic properties of the flow field that are suspected to play an important role in the onset of the instability (see (1.1)) if the instability is of the ‘curved streamlines’ type. Figure 10 illustrates the return to symmetry induced by the cylinder by performing increasing blockage ratio simulations at a constant value of $Wi=0.7$ with $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ .

In order to better characterize the magnitude of the observed instability in the numerical simulations, we define an asymmetry parameter as the ratio of the $x$ component of velocity to the mean value of velocity in the inlet arms, i.e. $AP=(\tilde{U}_{x})/\tilde{U}_{B}$ , at the $x=0,y=0.5$ location. Under this definition, if the flow retains its symmetric distribution, the $x$ component of velocity at this location will be equal to zero, while once the instability is triggered, due to the asymmetric distribution of the resulting flow, $\tilde{U}_{x}$ exhibits a non-zero value.

Figure 11. Effect of Weissenberg number on the symmetry-breaking instability for $\unicode[STIX]{x1D6FC}=0.02$ , $\unicode[STIX]{x1D6FD}=1/9$ and different blockage ratios $\unicode[STIX]{x1D6F7}$ .

Figure 12. Effect of blockage ratio parameter on the symmetry-breaking instability for $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ .

Previous studies (Poole et al. 2007a , b ) showed that by increasing the Weissenberg number, the symmetry breaking exhibits a supercritical growth close to the bifurcation point. In figure 11, the variation of the asymmetry parameter is plotted versus the Weissenberg number for different values of blockage ratios, which all exhibit a similar supercritical characteristic behaviour for the instability as was observed previously in the standard cross-slot (Poole et al. 2007b ). In figure 12, we show the effect of the blockage ratio for different constant values of the Weissenberg number with $\unicode[STIX]{x1D6FD}=1/9$ and $\unicode[STIX]{x1D6FC}=0.02$ . These results show a stabilizing effect of the addition of the cylinder which is characterized by a supercritical behaviour close to the bifurcation point $AP^{2}=a(1-\unicode[STIX]{x1D6F7})+b$ , where values of $a$ and $b$ are shown in figure 12.

Previous experimental studies with the standard cross-slot geometry (Arratia et al. 2006; Pathak & Hudson 2006; Haward et al. 2012a ; Sousa, Pinho & Alves 2018) have shown that by increasing the Weissenberg number to higher values, one can potentially trigger a second time-dependent instability. A numerical study of shear-thinning sPTT fluids conducted by Cruz et al. (2016) has also shown that the critical values of Weissenberg number for the onset on instability for both the steady symmetry-breaking and the time-dependent instabilities are a function of the cross-section aspect ratio ( $AR=\text{height}\text{/}\text{width}$ ) and the shear-thinning properties of the fluid (i.e. the $\unicode[STIX]{x1D6FC}$ parameter in the sPTT model). In Cruz et al. (2016), it was shown that once the aspect ratio is sufficiently small, depending on the value of the $\unicode[STIX]{x1D6FC}$ parameter, the flow distribution can potentially switch directly from the steady symmetry state to the time-dependent instability. In order to test our hypothesis that a cylinder can delay the steady symmetry-breaking instability, we extend our analysis from 2-D numerical simulations to a series of experiments using the flow-visualization approach introduced in § 3. For this purpose, we use two cross-slot geometries with cross-section aspect ratio of two, $AR=2$ , one standard cross-slot geometry and one cross-slot geometry with a cylinder at the geometrical centre with $\unicode[STIX]{x1D6F7}=0.55$ . This geometry is selected to best approximate a 2-D channel and observe the steady asymmetry in the standard geometry based on the results of Cruz et al. (2016).

Figure 13. Visualization of the flow pattern and schematic representation of characteristic lengths used in the definition of the asymmetry parameter in our experiment for the standard cross-slot geometry ( $\unicode[STIX]{x1D6F7}=0$ ) at $Re=6\times 10^{-6}$ for solutions of (a) 70 : 30 glycerine–water Newtonian fluid and (b) 190 ppm PAA in 70 : 30 glycerine–water viscoelastic fluid at $\bar{Wi}=0.8$ . Fluid with the fluorescent dye is injected from right inlet arm whereas the fluid without dye is injected from left.

A florescent dyed fluid is injected from the right inlet while undyed fluid is pumped in from the left inlet, both at the same flow rate, to visualize the flow distribution in the studied geometries. By increasing the flow rate, we are able to change the Weissenberg number from 0.1 up to 200. Our results suggest that beyond a critical value of the Weissenberg number, $\bar{Wi}_{cr}\approx 0.46$ , the symmetry of flow is broken in agreement with previous studies (Arratia et al. 2006; Sousa et al. 2018). In the standard cross-slot geometry, figure 13 shows the flow distribution in symmetric Newtonian flows and a steady asymmetric flow for 190 ppm PAA in 70 : 30 glycerine–water viscoelastic fluid at $\bar{Wi}=0.8$ and $Re=6\times 10^{-6}$ . To better investigate the symmetry-breaking instability we define an experimental asymmetry parameter as $AP_{e}=((\tilde{H}11-\tilde{H}12)+(\tilde{H}22-\tilde{H}21))/2\tilde{W}$ . A schematic definition of $\tilde{H}11$ , $\tilde{H}12$ , $\tilde{H}21$ and $\tilde{H}22$ is illustrated in figure 13(a).

Figure 14. Variation of the asymmetry parameter against the Weissenberg number for 190 ppm PAA in 70 : 30 glycerine–water viscoelastic fluid in standard cross-slot geometry.

Figure 15. Flow visualization at different instants of the oscillation cycle at the centre of the standard cross-slot geometry at $\bar{Wi}=277$ . Fluid with the fluorescent dye is injected from right inlet arm whereas the fluid without dye is injected from left.

Figure 16. Visualization of the flow pattern in cross-slot geometry with $\unicode[STIX]{x1D6F7}=0.55$ at $Re=2\times 10^{-5}$ for solutions of (a) 70 : 30 glycerine–water Newtonian fluid and (b) 190 ppm PAA in 70 : 30 glycerine–water viscoelastic fluid at $\bar{Wi}=2.3$ . Fluid with the fluorescent dye is injected from right inlet arm whereas the fluid without dye is injected from left.

Figure 17. Flow visualization at different instants of the oscillation cycle at the centre of the cross-slot geometry with $\unicode[STIX]{x1D6F7}=0.55$ and $\bar{Wi}=18.5$ . Fluid with the fluorescent dye is injected from right inlet arm whereas the fluid without dye is injected from left.

In figure 14, we show that the symmetry-breaking instability in the standard cross-slot geometry which, near the bifurcation point, follows a supercritical growth was also observed in our 2-D numerical simulations. As the Weissenberg number is increased to higher values, a time-dependent symmetry-breaking instability develops at $\bar{Wi}\approx 92$ which oscillates between two states (see figure 15) with a period of approximately $5\bar{\unicode[STIX]{x1D706}}~(\text{s})$ . Unfortunately, our experimental protocol was insufficient to characterize this time dependency more quantitatively.

Following the benchmark of our experimental protocol, we replace the standard cross-slot geometry with a modified geometry with $\unicode[STIX]{x1D6F7}=0.55$ . In figure 16, the symmetric distributions of the flow field for both the Newtonian fluid and 190 ppm PAA in 70 : 30 glycerine–water viscoelastic fluid at $\bar{Wi}=2.3$ and $Re=2\times 10^{-5}$ are presented. Here, unlike the standard cross-slot geometry where the purely elastic symmetry-breaking instability was triggered at $\bar{Wi}_{cr}=0.46$ , we were able to retain symmetric flow up to $\bar{Wi}_{cr}\approx 3.5$ , qualitatively supporting our numerical simulations regarding the stabilizing effect of the addition of a cylinder.

In the experiments, unlike the prediction from 2-D numerical simulations, we have noticed that by increasing the Weissenberg number to higher values, the flow field switches directly to an unsteady oscillatory instability (figure 17) with no evidence of the steady symmetry-breaking instability. This issue may potentially be related to the effect of finite cross-section aspect ratio in our experiments which was neglected in our 2-D simulations. Having shown experimentally that the cylinder does indeed stabilize the flow field significantly, we return to numerical simulations to probe additional effects and to see if we can quantify how the blockage ratio delays onset.

Figure 18 shows the ‘stability diagram’ for the modified geometry. As can be seen, above the black dashed line (the data are shown by red filled diamonds) are the regions in which the flow exhibits a steady asymmetric distribution. Our simulations suggest that as the size of the cylinder approaches zero, its modifying effects on the symmetry-breaking instability fade away and the critical value of the Weissenberg number for the onset of the instability approaches the same value as in the standard cross-slot geometry. One should note that a combination of conservation of mass (i.e. $(\unicode[STIX]{x2202}U_{r}/\unicode[STIX]{x2202}r)+(1/r)(\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703})=0$ ) with the no-slip boundary condition (i.e. velocity at the cylinder wall is equal to zero and due to the no-slip boundary condition $\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}|_{r=\unicode[STIX]{x1D6F7}}=0$ ) requires the strain rate at the stagnation point at the cylinder to become zero. So, by replacing the free stagnation point with pinned stagnation points, independent of the cylinder size, a non-zero value of strain rate at the stagnation point is replaced with a zero value once the no-slip boundary condition is applied (as confirmed in figure 19).

By changing the applied boundary condition at the cylinder surface at a constant blockage ratio, we attempt to analyse the onset sensitivity of the purely elastic instability with the magnitude of the strain rate in the vicinity of the cylinder. In figures 19 and 20, the strain rate distribution along the symmetry lines before the onset of instability is plotted using both slip and no-slip boundary conditions at the cylinder surface. The obtained results suggest that by replacing the no-slip boundary condition with the slip boundary condition, although the zero strain rate value at the stagnation point is replaced by a larger non-zero value (approximately 3.5 times larger than the maximum value of strain rate at the standard cross-slot geometry), the critical value of the Weissenberg number for which the instability is triggered stays almost constant at least in the small-aspect-ratio limit (i.e. $Wi_{cr}\approx 0.51$ for $\unicode[STIX]{x1D6F7}=0.05$ ).

Figure 18. Boundary between the symmetric and asymmetric flow for different $Wi$ and blockage ratio for $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ .

Figure 19. Effect of blockage ratio parameter on the flow strain rate before the start of instability for $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ along (a) horizontal line and (b) vertical line. (The critical values of $Wi$ for $\unicode[STIX]{x1D6F7}=0,0.25,0.40,0.50,0.55,0.60,0.65$ are $0.51,0.54,0.59,0.66,0.7,0.78,0.9$ , respectively.)

Figure 20. Effect of slip boundary condition at the cylinder on flow strain rate along (a) horizontal line and (b) vertical line for creeping Newtonian flows and viscoelastic fluids before the start of the purely elastic instability for $\unicode[STIX]{x1D6FC}=0.02$ , $\unicode[STIX]{x1D6FD}=1/9$ and $\unicode[STIX]{x1D6F7}=0.05$ .

Figure 21. Effect of (a) geometry modification and (b) slip boundary condition on the magnitude of the velocity field for creeping Newtonian fluid flows.

As shown in figure 21, the effects of the geometry modification and change of the boundary condition for $\unicode[STIX]{x1D6F7}=0.10$ are essentially local and mainly influence the flow distribution near the central region. Our results presented in figures 18–21 show that in the limit where the blockage ratio tends to zero ( $\unicode[STIX]{x1D6F7}\rightarrow 0$ ), even by applying fundamental changes on the kinematic characteristics of the flow field in the vicinity of the stagnation point (changing both the strain rate value and the curvature of streamlines), the critical value of the Weissenberg number for which the instability is triggered tends towards its critical value for the standard cross-slot geometry (i.e. $Wi_{cr}=0.51$ ). These data strongly suggest the symmetry-breaking instability in the cross-slot geometry is triggered at a location away from the geometric centre and closer to the corner regions as it is essentially independent of the kinematic properties of the flow near the stagnation point.

Figure 22. Illustrative contours of $M$ parameter for $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ .

Figure 23. Variation of critical values of Weissenberg number against the blockage ratio parameter for $\unicode[STIX]{x1D6FC}=0.02$ and $\unicode[STIX]{x1D6FD}=1/9$ (analytical fit: $1/Wi_{cr}=-1.9578\unicode[STIX]{x1D6F7}^{2}+1.9858$ ).

From the results presented we conclude that the high shear rate along with the high streamline curvature near the four cross-slot corners is more likely to be responsible for the transition to an asymmetric flow pattern and we therefore apply the ideas developed by McKinley et al. (1996) to investigate this possibility. Following the numerical procedure explained by Cruz et al. (2016), we have plotted the local distribution of the $M$ parameter (1.1) in figure 22. Local illustration of the $M$ parameter can be an effective way to observe the location of instability-driving regions within the flow field. As expected from previous results, the maximum value of the $M$ parameter appears near the corners which is related to the high curvature of streamlines and shear rate in this region (Hinch 1993).

One of the interesting properties of the $M$ parameter is related to its ability to scale purely elastic instabilities with respect to both rheological and geometrical properties of the problem (McKinley et al. 1996; Pakdel & McKinley 1998; Alves & Poole 2007; Zilz et al. 2012). Here, we use the approach suggested by McKinley et al. (1996) to model the effect of blockage ratio. An initial analysis suggests that neither $\tilde{D}$ nor $\tilde{W}$ provides a good approximation of the characteristic curvature of the streamlines in the cross-slot cylinder geometry, and in fact this characteristic length scale is influenced by both of these parameters. In appendix A, using an approximate analytical solution for creeping Newtonian flows, we show that the streamline curvature should scale as

(4.2) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{1}{\tilde{\Re }}}={\displaystyle \frac{1}{\tilde{W}}}(a+b\unicode[STIX]{x1D6F7}^{2}), & & \displaystyle\end{eqnarray}$$

where $a$ and $b$ are undetermined constants. The results presented in figure 22 clearly show that the highest value of the $M$ parameter appears in the vicinity of the corner’s shear-dominated regions. Assuming a steady-state purely shear flow for an Oldroyd-B fluid (to which the sPTT model approaches in the $\unicode[STIX]{x1D6FC}\rightarrow 0$ limit), one may scale the normal-stress component in (1.1) as

(4.3) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D70F}}_{ss}=2\tilde{\unicode[STIX]{x1D706}}\tilde{\unicode[STIX]{x1D702}}_{p}\tilde{\dot{\unicode[STIX]{x1D6FE}}}^{2}. & & \displaystyle\end{eqnarray}$$

Using $\tilde{U}_{b}$ and $\tilde{\Re }$ as references for the velocity field and the length scale, the reference shear rate may be expressed as $\tilde{\dot{\unicode[STIX]{x1D6FE}}}=\tilde{U}_{b}/\tilde{\Re }$ . Substitution of these equations into the dimensionless $M$ criteria (1.1) results in the following condition for the onset criteria of purely elastic instability in the modified cross-slot geometry:

(4.4) $$\begin{eqnarray}\displaystyle M_{c}=\sqrt{{\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}\tilde{U}_{b}}{\tilde{W}}}(a+b\unicode[STIX]{x1D6F7}^{2})2(1-\unicode[STIX]{x1D6FD}){\displaystyle \frac{\tilde{\unicode[STIX]{x1D706}}\tilde{U}_{b}}{\tilde{W}}}(a+b\unicode[STIX]{x1D6F7}^{2})}, & & \displaystyle\end{eqnarray}$$

where $a$ , $b$ and $M_{c}$ are unknown constants. However, factoring out the quantity $(\tilde{\unicode[STIX]{x1D706}}\tilde{U}_{b}/\tilde{W})^{2}$ and with some rearrangement one can simplify the equation into

(4.5) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{1}{Wi_{cr}}}=\bar{a}+\bar{b}\unicode[STIX]{x1D6F7}^{2}, & & \displaystyle\end{eqnarray}$$

where $\bar{a}=a\sqrt{2(1-\unicode[STIX]{x1D6FD})}/M_{c}$ and $\bar{b}=b\sqrt{2(1-\unicode[STIX]{x1D6FD})}/M_{c}$ are now the unknown constants. Figure 23 shows the fit obtained based on (4.5) and our 2-D numerical simulations. From this prediction (4.5), one can show $\unicode[STIX]{x2202}Wi_{cr}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}=0$ once $\unicode[STIX]{x1D6F7}\rightarrow 0$ , suggesting that, in the limit $\unicode[STIX]{x1D6F7}\rightarrow 0$ , the modifying effect of the addition of a cylinder on the critical kinematic properties of the flow field and the onset criteria of instability fades away (figure 18 also illustrates this point). Bearing in mind that (4.5) is obtained using a pure-shear flow assumption (i.e. (4.3)) and scaling of streamline curvature is obtained using a Taylor expansion in the vicinity of corners, the surprisingly good agreement between the analytical prediction (4.5) and numerical simulations once again suggests that the instability is triggered due to the shear flow near the corners and not the elongational-dominated flow at the stagnation point itself. Finally, to test the sensitivity of our chosen model parameters, in appendix B we show how the extensibility parameter in the sPTT model modifies the critical conditions and in appendix C we show that the same results can be observed with the FENE-P model in certain parameter limits.

5 Conclusions

In this work a passive control mechanism is introduced to the cross-slot geometry by the addition of a cylinder at the geometric centre of the domain. We use a series of numerical simulations and an approximate analytical analysis, using the definition of the $M$ parameter introduced by McKinley et al. (1996), to scale and analyse the stabilizing effect of the proposed geometrical modification on the onset criteria for a purely elastic instability and compare these results to our experimental data. Our work suggests that, in the limit $\unicode[STIX]{x1D6F7}\rightarrow 0$ , replacing the finite strain rate ‘free’ stagnation point flow with zero strain rate ‘pinned’ stagnation points at the cylinder surface does not affect the onset for the purely elastic instability. We extend our simulations by applying a complete slip boundary condition at the cylinder thereby changing the maximum strain rate from zero to approximately 3.5 times the strain rate in the standard cross-slot geometry, to show that the kinematic properties of the flow distribution around the stagnation point do not play any significant role in the onset of the purely elastic symmetry-breaking instability and that the critical values of the Weissenberg number for which the instability is triggered, in both cases, tend toward the critical value for the standard cross-slot geometry (i.e. $Wi\approx 0.51$ ). Finally, by plotting the local distribution of the $M$ parameter, we show that the location of instability-driving regions appears in the vicinity of the corners which can be attributed to the high deformation rate and strong streamline curvature in this region. We also show how the effect of blockage ratio scales with the onset criteria and thus can be successfully used to predict onset conditions.