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

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.


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), expansions (Pinho et al. 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 (Alves et al. 2001;McKinley 2002;Bisgaard 1983;Walters & Tanner 1992). 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 (Zilz et al. 2014;Poole et al. 2007a). During the last few decades, significant attention has been placed on such flows to characterise the physics behind these type of instabilities especially in shear flows (McKinley et al. 1996;Larson et al. 1990), however understanding of this type of elastic instability in extensional-dominated flows is not as advanced as for shear-dominated flows (Haward et al. 2016). The cross-slot geometry is a common flow geometry that has been utilized for generating controllable planar elongational flowfields and to study the stretching dynamics of polymers (Haward et al. 2016(Haward et al. , 2012b. This geometry consists of four bisecting rectangular channels with two sets of opposite inlets and outlets. These opposing inlets and outlets produce a flowfield 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 realised. To obtain an ideal planar elongational flowfield, Haward et al. (2012b) suggested to optimise the standard shape of the cross-slot geometry using a numerical approach and termed the resulting shape the Optimized Shape Cross-slot Elongation Rheometry geometry ("OSCER"). 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 (McKinley 2002;Bisgaard 1983;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 et al. 2005) and free stagnation point flows (Bisgaard 1983). Soulages et al. (2009) studied two different geometries to investigate the kinematic differences between a pinned stagnation point flow at the 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 3D 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 in between the corner and the stagnation point experience both significant shear flow near these re-entrant corners (Moffatt 1964;Dean & Montagnon 1949;Hinch 1993;Davies & Devlin 1993) and elongational-dominated flow near the stagnation point (Haward et al. 2016;Öztekin et al. 1997). The combination of this complex deformation with the non-linear 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 characterised. Supporting numerical simulations for this phenomenon were presented by Poole et al. (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 crossslot geometry this is debated (Wilson 2012) and something we will directly address in the current manuscript. 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: whereŨ is the magnitude of the local velocity,˜ is the local radius of curvature of a streamline,τ ss is the normal stress in the streamwise direction,η 0 is the zero shear rate viscosity of the fluid andγ is the magnitude of the shear rate. Throughout this paper we will indicate dimensional variables using a tilde. In equation (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 equation (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 in 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 the upper-convected Maxwell (UCM) and simplified Phan-Thien Tanner (sPTT) fluids. Previous studies for flow around sharp re-entrant corners both for Newtonian (Dean & Montagnon 1949;Moffatt 1964), and UCM fluids (Davies & Devlin 1993;Hinch 1993) 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 flowfield 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 if 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 et al. 2017).
In the current work, using a series of numerical simulations supported by experimental and analytical tools, we will 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 section 2 the governing equations and the employed numerical procedure is reviewed. In section 3, the experimental set-up and the employed protocol is explained. Following this, in section 4, a discussion on the obtained results in this study is presented followed by the conclusions.

Numerical procedure
The governing equations for this problem are conservation of mass, assuming incompressibility, and momentum: whereρ is the density of the fluid,ũ is the velocity vector,η s is the solvent viscosity andτ is the extra-stress tensor containing the polymeric contribution to the stress. In this work two different constitutive equations, namely the simplified Phan-Thien and Tanner (sPTT) (Phan-Thien & Tanner 1977) and the finitely extensible non-linear elastic (FENE-P) (Bird et al. 1980) models, have been used to study the effect of different constitutive equations.
The extra-stress tensor using the sPTT model may be calculated as follows: whereλ is the relaxation time of the fluid andη p is the polymeric contribution to the viscosity. The upper-convective derivative of the extra-stress tensor, ∇ τ , is defined as: where the material derivative of an arbitrary matrixÃ is defined as D Dt (Ã) = ( ∂Ã ∂t + u.∇Ã). The f 1 function for the linear-sPTT model is defined as follows: where α is the the extensibility parameter and in the limiting case that α = 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 et al. 1987).
One other constitutive equation which has frequently been used to simulate the shearthinning behaviour of polymeric materials is the FENE-P (Bird et al. 1980) model which may be presented in the following form: where a = L 2 L 2 −3 and L 2 is also called the extensibility parameter. The f 2 function for the FENE-P model is defined as: (2.7) 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 globalx,ỹ andz coordinate system, the logarithm 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 the works conducted by Afonso et al. (2012) and Fattal & Kupferman (2004).

Problem definition
A schematic of the geometry and employed coordinate systems are shown in Figure  1. The length of inlet/outlet arms are set to be 15 times of 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Ũ B at the inlets, and Neumann boundary condition at the outlets were applied. At the walls, the 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.

Non-dimensionalization
It is convenient to use dimensionless parameters in this problem. The relationships between dimensional and dimensionless parameters are: werex,ỹ andr are the variables related to the corresponding rectangular and polar coordinate systems,W is the width of the channel,D is the diameter of the cylinder, Φ is the blockage ratio parameter which may change between zero (the standard cross-slot geometry with no cylinder) to √ 2 (fully blocked cross-slot geometry) values,Ũ is the velocity vector,Ũ 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, W i is the Weissenberg number, β is the solvent-to-total viscosity ratio andη t is the total viscosity (i.e.η t =η s +η p ), N 1 is the non-dimensional first normal-stress difference,˜ is the strain rate (either˜ xx or˜ yy ) andΩ is the vorticity tensor.

Mesh dependency study
In this section, a number of representative results analysing the effect of mesh on the 2D 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 sixteen 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 Φ = 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 (GR) 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.
In Figure 3, the effect of mesh refinement on the flow distribution of creeping Newtonian flow are presented (i.e. Re = 0.01). In this figure, the distribution of velocity magnitude at the entrance region of the outlet arms, are 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 Φ = 0.10 in which meshes similar in density to M2 were used. For the standard cross-slot geometry a mesh similar to the work by

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 (m) certified accuracy is 0.35% at the lowest free pulsation-free delivering rate, which set our lowest flow rate.
For the purpose of this study, two geometries including one standard cross-slot device and one modified cross-slot geometry with Φ = 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 2D channel, the channel height is selected to be twice of its width,H = 2000 ± 10[µm] andW = 1000 ± 10[µm], respectively. The device was encased in polyoxymethylene, an insulating material also known as ACETAL. The cross-section dimensions of the cross-slot flow device were quantified using a Nikon EPIPHOT TME inverted microscope with 100 times magnification, 470 pixels = 1000 [µm]. The combined length of the channel inlet and outlet is 80W = 80 ± 0.1[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 visualised.
Rhodamine-B (ACROS Organics) was chosen as a fluorescent dye to capture the flow patterns in the cross-slot micro-geometry (Ross et al. 2001;Huang et al. 2013Huang et al. , 2014. The  . (a) A schematic illustrating the experimental apparatus of a microfluidic cross-slot device allowed 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 [notes: 1-not to scale; 2-camera FOV to be captured (field of view): 4mm × 4mm]. (b) The prototype microfluidic device rigs for I: standard cross-slot geometry and II: the cross-slot with cylinder geometry. The channels were micro-machined in brass and encased in polyoxymethylene. (c) Photograph illustrating the experimental rig set-up assembled.
to excite this fluorescent Rhodamine-B dye with illumination and a CCD (charge couple device) camera.

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 percent glycerine (by weight). The density ( water (by weight) and 190 ppm polyacrylamide (PAA, M w = 18 × 10 6 [g.mol −1 ], Sigma-Aldrich). 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 hours, the result is two clear, transparent, and colourless solutions. The PAA and glycerine solutions are then combined in a unique container, and the resultant polymer solution is thus left to mix gently in a magnetic stirrer (24 hours). All rheology measurements were conducted atT = 20 ± 0.1[ • C]. The working fluid characterisation 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 α = 0.63,η 0 = 3.36[P a.s] andη ∞ = 35.5[mP a.s]. The relaxation time of the fluid were measured using small amplitude oscillatory shear (SAOS). Analysing the storage modulus (G ) and the loss modulus (G ) using a Maxwell model is a standard approach to estimate 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:

asλ
As for the shear viscosity, the Anton Paar MCR302 is used to measure the SAOS properties of the working fluid. Both the amplitude and frequency sweep tests are carried out at T = 20[ • C], using a 60[mm] and 1 • cone. Figure 6 shows the results of the frequency sweep test in the range of 0.01 to 4 [rad/s] at 10% strain amplitude (which we confirmed was in the linear limit). Our experimental results suggest an average relaxation time of λ = 45.05[s].

Experimental protocol
Experiments are conducted over a range of 0.1 <W i < 200 by programming the syringe pumps to perform ramps inW i with small step increases in volumetric flow rate (0.001 <Q < 0.2[nl/s]). Here,W i is defined based on the average relaxation time of the fluid obtained from equation (3.3). The camera is synchronised with the laser at a repetition rate of 8.875 [Hz]. All images shown in this work were captured from the top view using an 8× microscope objective (Leica Microsystems GmbH), and the camera FOV (field of view), 4mm x 4mm which thoroughly covers the central region of the device containing the field of interest. The flow visualisation 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 but two identical syringes. Fluid motion starts from rest and at least fifty 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 assure 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.

Results and Discussion
In this section, results related to the influence of the geometric modification on the flowfield and on the associated steady symmetry-breaking instability in the cross-slot geometry are presented. We use a series of 2D numerical simulations, experimental results and a supporting approximate analytical solution to determine the M parameter (equation 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. 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 visualise different types of flow is the flow-type parameter ξ (Lee et al. 2007):  Figure 7 show the effect of elastic stress on the streamline distribution superimposed on the flow-type parameter for small values of W i number before the onset of the symmetry-breaking purely-elastic instability. In all figures, the inlets are located on the left and right sides while the outlets are on 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 flowfield in this area. As shown in Figure 7, in the creeping flow of Newtonian fluids, this extensionaldominated 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 (Ω xy ) as shown in the right hand column of Figure 7.
By increasing the Weissenberg number up to a certain limit (here, W i = 0.17 for β = 1/9 and α = 0.02), the velocity distribution in the fully-developed inlet arms exhibits a flatter distribution in comparison to the Newtonian flowfield, due to shearthinning, and consequently the elongational dominated flow expands to a wider region. As shown in Figure 8, by increasing the value of 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 = ±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 flowfield, even prior to any purely-elastic instability, is significant and changes the flow-type downstream of the stagnation point in a quite fundamental way.
Previous experimental studies (Arratia et al. 2006;Gardner et al. 1982) conducted in the standard cross-slot geometry have shown that beyond a critical value of the W i number, 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 symmetry-breaking instability is presented in Figure 9 (for a nominal case of Φ = 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. By increasing the blockage ratio, the average velocity passing through the central region of the cross-slot geometry (i.e. throughL 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 to control the kinematic properties of the flowfield that are suspected to play an important role in the onset of the instability (see equation 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 W i = 0.7 with α = 0.02, β = 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 , 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,Ũ x exhibits a non-zero value.
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 β = 1/9, α = 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 − Φ) + b where values of a and b are shown in Figure 12. Previous experimental studies in the standard cross-slot geometry (Arratia et al. 2006;Pathak & Hudson 2006;Haward et al. 2012a;Sousa et al. 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 = height width ) and the shearthinning properties of the fluid (i.e. the α 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 α parameter, the flow distribution can potentially switch from the steady symmetry state to the time-dependent instability directly. In order to test our hypothesis that a cylinder can delay the steady symmetry-breaking instability, we extend our analysis from 2D numerical simulations to a series of experiments using the flow-visualization approach introduced in Section 3. For this purpose, we use two cross-slot geometries with the crosssection aspect ratio of two, AR = 2, one standard cross-slot geometry and one cross-slot geometry with a cylinder at the geometrical centre with Φ = 0.55. This geometry is selected to best approximate a 2D channel and observe the steady asymmetry in the standard geometry based on the results of Cruz et al. (2016).
A florescence 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 visualise 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 suggests that beyond a critical value of the Weis- senberg number,W i cr ≈ 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 asymmetry flow for 190 ppm polyacrylamide (PAA) in 70:30 glycerine/water viscoelastic fluids at W i = 0.8 and Re = 6e − 6. To better investigate the symmetry-breaking instability we define an experimental asymmetry parameter as AP e = (H11−H12)+(H22−H21)
In Figure 14, we show that the symmetry-breaking instability in the standard crossslot geometry which, near the bifurcation point, follows a supercritical growth as was also observed in our 2D numerical simulations. As the Weissenberg number is increased to higher values, a time-dependent symmetry-breaking instability develops atW i ≈ 92 which oscillates between two states (see Figure 15) with a period of, approximately, 5λ[s]. Unfortunately, our experimental protocol was insufficient to characterise this time dependency more quantitatively. Following the benchmark of our experimental protocol, we replace the standard crossslot geometry with a modified geometry with Φ = 0.55. In Figure 16, the symmetric distribution of the flowfield for both the Newtonian fluid and 190 ppm polyacrylamide (PAA) in 70:30 glycerine/water viscoelastic fluids atW i = 2.3 and Re = 2e − 5 are presented. Here, unlike the standard cross-slot geometry where the purely-elastic symmetry-breaking instability was triggered atW i cr = 0.46, we were able to retain symmetric flow up toW i cr ≈ 3.5, qualitatively supporting our numerical simulations regarding the stabilizing effect of the addition of a cylinder.
In the experiments, unlike the prediction from 2D numerical simulations, we have noticed that by increasing the Weissenberg number to higher values, the flowfield switches to an unsteady oscillatory instability ( Figure 17) directly 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 2D simulations. Having shown experimentally that the cylinder does indeed stabilise the flowfield 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 to zero, its modifying effects on the symmetry-breaking instability fades away and the critical value of the Weissenberg number for the onset of the instability approaches to the same value as in the standard cross-slot geometry. One should note that a combination of conservation mass (i.e. ∂Ur ∂r + 1 r ∂U θ ∂θ = 0) with the noslip boundary condition (i.e. velocity at the cylinder wall is equal to zero and due to the no-slip boundary condition ∂U θ ∂θ | r=Φ = 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 at 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 suggests 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 in which the instability is triggered stays almost constant at least in the small aspect ratio limit (i.e. W i cr ≈ 0.51 for Φ = 0.05).
As shown in Figure 21, the effects of the geometry modification and change of the boundary condition for Φ = 0.10, is essentially local and mainly influences 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 (Φ → 0), even by applying fundamental changes on the kinematic characteristics of the flowfield 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 in which the instability is triggered tends towards its critical value for the standard cross-slot geometry (i.e. W i 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.
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 (equation 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 flowfield. 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 rhelogiocal and geometrical properties of the problem (Pakdel & McKinley 1998;McKinley et al. 1996;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 neitherD norW provide 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: 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 (of which the sPTT model approaches in the α → 0 limit), one may scale the normal-stress component in equation 1.1 as: UsingŨ b and˜ as references for the velocity field and the length scale, the reference shear rate may be expressed asγ =Ũ b . Substitution of these equations into the  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) dimensionless M criteria (equation 1.1), results in the following condition for the onset criteria of purely-elastic instability in the modified cross-slot geometry: where a, b and M c are unknown constants. However, factoring out the quantity (λŨ b W ) 2 and some rearrangement one can simplify the equation into: 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 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.

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 that Φ → 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 of 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 in which the instability is triggered, in both cases, tend toward its critical value for the standard cross-slot geometry (i.e. W i ≈ 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. In this section, an approximate analytical solution for the flow distribution of a creeping Newtonian fluid along a diagonal line between a corner and a stagnation point is obtained (a schematic of the problem is depicted in Figure 1). Moving to a cylindrical coordinate system, the components of the velocity vector in terms of the stream function are: For steady Newtonian flows, the conservation of momentum equation using the definition of stream function may be presented as: Separation of the stream function into a function of a radial line f (r) and a function of the polar angle which is periodic with the wave number m allows us to obtain the general form of the analytical solution for equation A 2 as: The analytical solution obtained by Moffatt (1964) for the creeping flow of Newtonian fluids around a sharp corner implies that the stream function must have the following form: 28ψ = Ar 1.5445 (cos(0.3415π)cos(1.5445θ) − cos(1.158π)cos(0.4555θ)).

(A 4)
Using the general form of the stream function solution (equation A 3), we aim to find an approximate analytical expression for the stream function which as r → 0 the solution asymptotes to the exact solution presented by Moffatt (1964) (i.e. equation A 4). One should note that 1 < m < 3 , which implies that the velocity is not singular atr = 0 and the velocity gradient has a bounded value far from the corner. Equation A 4 shows that asr → 0 the stream function must have wave numbers (the m parameter in equation A 3) of 0.4555 and 1.5445. Using the suggested values for m in equation A 3, the stream function is now: ψ = Ar 1.5445 (cos(0.3415π)cos(1.5445θ)−cos(1.158π)cos(0.4555θ))+Br 2.4555 cos(0.4555θ).
(A 5) considering the following constraints at θ = 0: The unknown constants can then be calculated as: A comparison between this approximate analytical solution and the numerical result is shown in Figure 24 where a good agreement can be observed. Knowing the functional form of the velocity (equation A 5), the curvature of the streamline may be calculated as (Haward et al. 2016;Cruz et al. 2016): where Dũ Dt is the material derivative of the velocity vector. Using the definition of the stream function, one can calculate 1/ as follows: 1 = 1.538(L) 0.911 − 3.811r 0.911 (4.163(L) 0.911 − 4.163r 0.911 )r .
(A 11) One should notice that the present analytical solution is only an approximate solution for creeping Newtonian fluid flow between the corner of the cross-slot geometry and the stagnation point and so the obtained stream function does not satisfy the boundary conditions at the walls and does not present a correct θ dependency away from θ = 0.
To include the effect of the cylinder at the center of the cross-slot geometry, we consider the diagonal gap lengthL to be simply: Using equation A 12 and some simplification with a Taylor expansion around the corner of the cross-slot geometry, the 1/˜ term can be expressed as: Equation A 13 suggests a second order dependency of the streamline curvature with the blockage ratio parameter. Here we assume that, prior to instability, viscoelasticity does not significantly affect the scaling correlation between the streamline curvature and the blockage ratio and will proceed to scale the M parameter with the blockage ratio parameter based on equation A 13. In the end we will check the accuracy of this hypothesis by plotting the obtained approximate analytical expression with our numerical results.
For an Oldroyd-B fluid, assuming steady-state simple shear, the second term on the right hand side of equation 1.1 can be simplified to: ConsideringŨ B as our reference velocity, and using the derived˜ parameter as our reference length (i.e. equation 4.2), after some simplifications, one can state: so that a plot of the reciprocal of the critical condition (1/W i c ) against Φ 2 should be linear with intercept a * and slope b * (Figure 23 shows this plot confirming the approximate analysis). (1 +λ (T r(τ )) η p L 2 )τ +λ ∇ τ =η p (∇ũ + ∇ũ T ).

(C 2)
This way of expressing the FENE-P model suggests that in "viscometric" flows, assuming L 2 >> 1, both the FENE-P and sPTT models are identical when α = 1 L 2 (Cf Equation C 2 and equations 2.3 and 2.5). In Figure 27, this point is illustrated by plotting viscometric functions of both the sPTT and FENE-P models with α = 2e − 4 and L 2 = 5000, respectively, with β = 1/9. Our simulations with L 2 = 5000 and α = 2e − 4 shown in Figure 28 suggests that once the viscometric functions are matched the instability is triggered at almost the same critical value of the Weissenberg number which may, at first glance, suggest that the critical condition in which the instability is triggered is related to the steady-state viscometric properties of the model. However, one should note that very large values of L 2 in the FENE-P and very small values of α in the sPTT model may potentially conceal the difference between these two models since in the limit L 2 → ∞ or α → 0 both the FENE-P and sPTT models will reduce to the Oldroyd-B model (Bird et al. 1987).
In another attempt, to better investigate the effect of rheological properties of the fluids, we try to fit the viscometric functions using a smaller value of L 2 and a larger value of the α parameter in these two models. Figure 29, shows the fit obtained between viscometric function of SPTT model with α = 0.02 and FENE-P model with L 2 = 58 at a fixed value of β = 0.2. Interestingly, as shown in Figure 30, although the steady-state viscometric functions between the two constitutive equations are in good agreement, the critical values of the Weissenberg number for the onset of the instability are significantly different in these two models. This result provides further evidence that the kinematic properties of the flow triggering the instability are, most probably, related to the nonhomogeneous aspects of this complex flow and not solely related to the stagnation point.