Polymer solution flow transitions and scaling laws for changing contraction ratios in planar constriction microchannels

Abstract

Pores scale flows through contractions and expansions are relevant in geoengineering, microfluidics and material processing etc. These flows experience shearing and extensional kinematics near constrictions, where polymer solutions may demonstrate instabilities that arise from the fluid's nonlinear rheological characteristics even in creeping flows. The relative effect of shearing and extension can be controlled by the flow geometry. Following our earlier reports on the constriction length (M. K. Raihan et al., Soft Matter, 2021, 17, 9198–9209) and depth (M. K. Raihan et al., Soft Matter, 2022, 18, 7427–7440), we investigate here the flow responses to changing constriction width and in turn contraction ratio, CR, of the main channel width to the constriction width in planar constriction microchannels. We test water and three polymer solutions including shear thinning xanthan gum, viscoelastic polyethylene oxide (PEO), and shear thinning/viscoelastic polyacrylamide solutions. Overall, the contraction and expansion flows in all tested fluids demonstrate destabilization with increasing CR except for the PEO solution, where the threshold Reynolds number for the onset of contraction flow instability first increases and then decreases. Such nonmonotonic CR dependence is also observed from the vortex length in the contraction PEO flow. In contrast, the vortex length for every other case has a fixed-order (either zero or a positive number based on fluid rheology) dependence on CR. The insights obtained here will benefit the designing of lab-on-a-chip devices as well as the harnessing of pore-scale flows for enhanced mixing, material recovery and sequestration purposes.

---

1 Introduction

Contraction and expansion microchannels are often encountered in porous media flows, lab-on-a-chip assays, inkjet printers, and material extrusion processes etc.^1–6 Near the constriction of such channels, the flow field undergoes considerable shearing and extensional kinematics. Under these conditions, polymer solutions can induce instabilities starting in the vicinity of any re-entrant or salient corners depending on the fluid rheology even in creeping flows.^7–12 This arises from the polymer chains transitioning from shear distorted to stretched conformations while passing through the constriction. This chain conformational dynamics manifests as flow instabilities based on how favorable the transitions are.^13–20 The ratio of shearing to extensional components in the flow field plays a key role in the polymer chain transitions. It can be directly varied by changing the contraction ratio, CR, of the main channel width to the constriction width.^4,5,8,9 Consequently, a wide number of studies in the non-Newtonian fluid mechanics community have investigated CR effects. These studies have typically considered contraction and expansion domains on this topic can be classified into two geometric groups, i.e., planar and axisymmetric constriction channels. The latter can incur more complexities due to the three-dimensionality of the constriction. Nevertheless, both these domains have been tackled in synergy over the past fifty years.

The earliest report on CR effects from Nguyen and Boger²¹ considered axisymmetric contraction channels, where increasing the CR was observed to enhance the instabilities in viscoelastic flows with pronounced feature asymmetries. A few years later, Evans and Walters²² noted the CR based secondary flow enhancements further in square–square and planar channels. Such CR dependent flow evolution patterns were demonstrated to vary also with the constriction type, where strong shear thinning was found to help maintain the stable symmetry of vortices. In a following investigation with a planar contraction channel, Evans and Walters²³ elucidated that changing the channel CR, polymer concentration, flow inertia, and the re-entrant corner rounding can result in vortex developments based at either the lips or the salient corners. Another experimental study from Chiba et al.²⁴ confirmed that the vortices in the planar contraction flow of shear thinning polyacrylamide (PAA) solutions start from the re-entrant corners and the onsetting flow rate decreases as CR increases. McKinley et al.²⁵ systematically reinvestigated the axisymmetric contraction flows with well characterized Boger fluids showing similar lip onsets and CR dependent amplifications of instabilities. They noted highly unstable asymmetric features forming only in relatively lower CR channels.

Emphasis was then extended to reproducing the flow events with simulations. The convenience of interpretation and computational costs have often motivated the numerical investigations to be performed in planar, if not simply two-dimensional, domains. Despite limitations in earlier results with the upper-convected Maxwell model,²⁶ Purnode and Crochet²⁷ reported successful qualitative agreements on the CR effects in viscoelastic contraction flows with the FENE-P model in a 2D finite element simulation. They investigated the influences of PAA concentration and re-entrant corner rounding under creeping flows that were experimentally studied by Evans and Walters.^22,23 With the advent of finite volume solver stabilization schemes, Alves et al.²⁸ carried out 2D simulations of viscoelastic flow through abrupt contractions with the Phan–Thien–Tanner (PTT) model for a range of CR and Deborah number, De (defined as the ratio of the fluid relaxation time to the characteristic time of the observation). They probed the purely viscoelastic features and found that the lip vortices scale with De while the salient corner vortices scale with De/CR. They also underlined the lip vortex streamlines and recirculation intensities being dominated by the downstream flow quantities.

Concurrently in the early 2000s, progress in soft lithography brought about the prospect of constructing resource efficient extensional micro-rheometers based on instability initiations and pressure drops. In that spirit, Oliveira et al.²⁹ studied experimentally the effect of extensional Hencky strain [defined as the natural logarithm of the ratio of the channel width to the minimum contraction width, i.e., ln(CR)] on the flow of water and a Boger solution of polyethylene oxide (PEO) in planar hyperbolic contraction-sudden expansion channels. They remarked that the pressure drop for the polymeric fluid is a strongly nonlinear function of Reynolds number, which poses difficulty in constructing the intended framework for extensional viscometry. The same group of authors later simplified the approach by considering only water flow features to experimentally elucidate the viscous effects along with numerical agreements.³⁰ Then in the same channel, Sousa et al.³¹ noted a monotonic destabilizing effect on the flow of aqueous xanthan gum (XG) solution, as a blood analog, with increasing CR. Another key takeaway of this study is that blood plasma must not be treated as a simple solution, especially in the micro-realm.

Based on the progress in understanding of planar domains, Oliveira et al.³² studied the CR effects in axisymmetric contraction channels using the Oldroyd-B and PTT constitutive models. They found similar dependencies of lip and salient corner features on De and De/CR, respectively, akin to the planar geometries.²⁸ In a following study, Sousa et al.³³ considered square–square contraction channels. Seven different regimes were identified for the Newtonian and Boger flows in the De−CR space with CR induced prominent changes to the flow occurring at higher De. Further considering a shear-thinning viscoelastic PAA solution, Sousa et al.³⁴ later charted four major regimes in their experiments. An interesting flow reversal phenomenon was noted in the vortices of polymer solutions contrasting the Newtonian case. Moreover, their numerical models with the PTT constitutive law, though for a limited range of De due to computational constraints, demonstrated the ability to capture the helical flow signatures of this domain, albeit qualitatively as the critical De for the instability onset was an order of magnitude lower than the experimental observation.

Recently Lanzaro and Yuan³⁵ investigated the CR effects on the flow of PAA solutions in planar abrupt constriction microchannels. They demonstrated that in addition to the Weissenberg and Renolds numbers, the viscoelastic flow phenomena are dependent on both the CR and aspect ratio of the flow geometry. Pérez-Camacho et al.³⁶ emphasized the combined effects of strong shear thinning and viscoelasticity in axisymmetric domains. The flow visualization revealed the size of corner vortices is related to the relative contributions from the first normal stress difference and extensional stress, independent of the nature of the fluid. In a following study, López-Aguilar et al.³⁷ used the so-called swanINNFM model to simulate the experiment from Pérez-Camacho et al.,³⁶ which was found capable of capturing the enhanced pressure drops in all CR channels. More recently, Hidema et al.³⁸ experimentally investigated the effects of polymer flexibility and entanglement on the elastic instability of sodium hyaluronate solutions in planar single contraction-expansion microchannels. Monotonic flow destabilization by increasing CR was observed without any major differences in feature developments. In a very recent study, Yin et al.³⁹ studied the flow regime transitions in PEO solutions with varying polymer concentrations through continuous abrupt contraction-expansion arrays with different CR. Their results indicate that a higher CR results in greater extensional rates, leading to faster polymer stretching, more intense polymer session, and quicker flow regime transitions.

Table 1 presents a summary of the experimental or numerical specifics in the above reviewed papers on polymer solution flows through constriction channels with changing CR. It is understood from these studies that the stability criteria and subsequent unstable feature developments with inertia in polymeric flows follow complex relationships with channel geometry and CR. In particular, the role of fluid rheology in the flow instabilities under considerable inertia for changing CR is not clear despite being of practical interest. Our group previously investigated how the geometric alterations of planar constriction, including length⁴⁰ and depth,⁴¹ affect the flow instabilities in shear thinning XG, viscoelastic PEO, and shear thinning/viscoelastic PAA solutions compared to Newtonian water. A systematic rheology-based knowledge with these fluids in planar constriction channels having different constriction widths would extend the current understanding of flow responses under these channel conditions which are, however, based on purely viscoelastic fluids. Therefore, we consider three constriction microchannels that vary in the constriction width and hence CR but similar in dimensions otherwise to experimentally investigate the polymer solution flow responses. We consider both creeping and moderately inertial flows yielding shear rate order of magnitude ranging between 10² s⁻¹ to 10⁵ s⁻¹ due to relevance in oil recovery, physiological flows, industrial processing, and lab-on-a-chip applications.^3–12 We summarize the flow regimes and self-similarities in flow instabilities and features observed due to CR alterations. It is expected that our findings in this work will be useful for lab-on-a-chip applications as well as understanding polymeric flows in porous media that pertain to widely ranging applications in fluid mixing, transport and reaction.

Table 1 Summary of previous studies on CR effects in polymer solution flows in axisymmetric (A), planar (P) and square (S) geometries. Note the Reynolds number, Re, and Weissenberg number, Wi, are defined in Section 2.2

Authors	Year	Exp./num.	Axisym./planar/square	CR	Fluid	Wi or De range	Re range	Comments
Nguyen and Boger^21	1979	E	A	4.09, 7.67, 12.3, 14.83	Glucose syrup, MCY41N, SWEETOL (Boger and shear-thinning)	0.13 < Wi < 0.30	10^−4 < Re < 10^−2	Increasing CR enhances instabilities in viscoelastic flows with pronounced feature asymmetries
Evans and Walters^22	1986	E	P, S	4, 16, 80	PAA (Boger and shear-thinning)	1.8 × 10^−6 < Wi < 9 × 10^−4	3 × 10^−5 < Re < 0.013	CR dependent instabilities vary with the constriction type and the reentrant corner rounding
Chiba et al.^24	1990	E	P	3.33, 10	PAA (shear-thinning)	No relaxation time	2.5 < Re < 133	Contraction vortices start from reentrant corners and CR monotonically enhances the process
McKinley et al.^25	1991	E	A	2, 3, 4, 5, 6, 8	PIB (Boger)	1.5 < De < 5.4	0.004 < Re < 0.137	An unstable asymmetric flow state occurs only at lower CR channels
Purnode and Crochet^27	1996	N	P	4, 16	FENE-P (shear-thinning)	0.033 < Wi < 145.2	3 × 10^−5 < Re < 6.37	Qualitative agreement with experiments of Evans and Walters^22 was achieved
Alves et al.^28	2004	N	P	4, 10, 20, 40, 100	PTT (shear-thinning)	1 < De < 250	Re ≪ 1	Lip vortices scale with De while corner vortices scale with De/CR
Oliveira et al.^29	2006	E	P	2.72, 7.39, 20.1, 40	PEO (Boger)	0.75 < De < 12.5	13.2 < Re < 145	Extra pressure drops vary highly nonlinearly with CR
Sousa et al.^31	2011	E	P	2.14, 6.42, 15.64	XG (shear-thinning)	59.8 < De < 732	3.47 < Re < 93	Blood plasma (modeled by XG) shows complex responses to CR changes in microchannels
Oliveira et al.^32	2007	N	A	2, 4, 10, 20, 40, 100	Oldroyd-B, PTT (Boger, shear-thinning)	0.01 < De < 1000	Re ≪ 1	Similar De and De/CR based scaling trends observed like Alves et al.^28
Sousa et al.^33	2009	E	S	2.4,4, 8, 12	PAA (Boger)	3 < De < 38	0.0148 < Re < 1.12	Square shape caused deviations in the vortex scaling trends seen in planar and axisymmetric channels for CR changes
Sousa et al.^34	2011	E	S	2.4, 4, 8, 12	PAA (shear-thinning)	0.0554 < De < 450	Re < 1	Viscoelasticity brings out an anomalous diverging flow upstream of the contraction
Lanzaro and Yuan^35	2011	E	P	4, 8, 16	PEO, PAA (Boger and shear-thinning)	1.4 < Wi < 131.7	0.02 < Re < 9.6	Flow responses to CR and aspect ratio of the channel also depend on the polydispersity of the samples
Pérez-Camacho et al.^36	2015	E, N	A	2, 4, 6, 8, 10	PAA (Boger), acrysol TT-935 (shear-thinning), FENE-CR	0.002 < Wi < 400	Re ≪ 1	Extra pressure drop responds monotonically to CR for Boger flows, but nonmonotonically for shear-thinning flows
Lopez-Aguilar et al.^37	2016	E	A	2, 4, 6, 8, 10	swanINNFM (Boger)	0.003 < De < 17.4	Re < 0.1	Extra pressure drops correlated well with vortices and extensional viscosity under different CR
Hidema et al.^38	2019	E	P	8, 16	Na-HA (Boger and shear-thinning)	1 < Wi < 1000	0.02 < Re < 20	Flows destabilized monotonically with CR
Yin et al.^39	2024	E	P	8, 16	PEO (Boger)	4 < Wi < 620	4.2 < Re < 40	CR is the most critical factor determining flow regime transitions

---

2 Experiment

2.1 Materials

Fig. 1(a) shows the top-view pictures of the contraction-expansion region of three single-constriction microchannel molds, which were fabricated with photo lithography using SU-8 photoresist and differed only in the width of the constriction. The microchannels were prepared with soft lithography based on polydimethylsiloxane and glass slides. They were each 500 μm wide and 2 cm long with a 250 μm long constriction in the middle. The width of the constriction ranged from 100 to 50 and 31 μm, yielding the contraction ratios, CR = 5, 10 and 16, respectively. All prepared channels were 50 μm deep everywhere.

Fig. 1 (a) Top-view pictures of the contraction-expansion region of constriction microchannel molds (important dimensions are highlighted) with the contraction ratio, CR = 5, 10 and 16, respectively; (b) experimentally measured (symbols) and theoretically fitted (lines, via the Carreau model) viscosity data of the prepared polymer solutions, where the vertical dashed line highlights the characteristic shear rate inside the constriction for the flow rate of 1 ml h⁻¹ in the CR = 10 constriction channel.

Three types of polymer solutions were prepared in DI water, which were reported to have distinct rheological properties: shear thinning 2000 ppm XG solution (molecular weight, M_w ≈ 2 MDa, Tokyo Chemical Industries), viscoelastic 1000 ppm PEO solution (M_w = 2 MDa, Sigma-Aldrich), shear thinning and viscoelastic 200 ppm PAA solution (M_w = 18 MDa, Polysciences).^42–44 Newtonian DI water was also tested as the control experiment. The shear viscosities of polymer solutions were measured using a cone-plate torsional rheometer (Anton Paar), which are presented in Fig. 1(b). The Carreau model was used to fit the viscosity data of shear thinning XG and PAA solutions,

where η is the dynamic viscosity, η_∞ is the infinite-shear-rate viscosity, η₀ is the zero-shear-rate viscosity, λ_c is a time constant expressed as the inverse of the critical shear rate when viscosity starts decreasing significantly, [small gamma, Greek, dot above] is the shear rate, and n is the power-law index. The viscosity of PEO solution varies negligibly with [small gamma, Greek, dot above] compared to the XG and PAA solutions and is assumed to be the average of the data points over the whole [small gamma, Greek, dot above] range. The relaxation times of the viscoelastic PEO and PAA solutions were extracted from previous studies^43,44 that reported the experimental data of the same polymer solutions as used here. The fitting parameters in the Carreau model, i.e., eqn (1), along with the relaxation times, λ, of the prepared fluids are summarized in Table 2.

Table 2 Summary of the rheological properties of the prepared fluids. The elasticity number, El (definition provided later in Section 2.2), was estimated at the flow rate of 10 ml h⁻¹ in the 10 : 1 constriction microchannel for all fluids

Solution	η 0 (mPa s)	η ∞ (mPa s)	λ C (s)	n	λ (ms)	El
a Rodd and coworkers.^43 b Poole and escudeier.^44
DI water	1	1	—	1	0	0
2000 ppm XG	1740	1.8	6.6	0.33	∼0	∼0
1000 ppm PEO	2.18	2.18	—	∼1	1.5^a	2.62
200 ppm PAA	4900	2.1	151	0.25	95^b	162

---

2.2 Methods

Spherical green-fluorescent particles of 1 μm diameter (Bangs Laboratories) were seeded in each of the prepared fluids at a concentration of 0.05% solid content. An infusion syringe pump (KD Scientific, Holliston, MA, USA) was used to drive the flow of particulate solutions through the three constriction microchannels. The top-view images of the tracer particles were recorded at the contraction-extraction region using an inverted fluorescent microscope (Nikon Eclipse TE2000U) with a CCD camera (Nikon DS-Qi1Mc) via a 10× objective lens. The exposure time of the camera was adjusted from 0.5 to 3 s for the captured snapshots, where the lower values were used to adapt to the faster movement of tracer particles. The Nikon imaging software (NIS-Elements AR 3.22) was used to acquire digital images and measure the size of vortices if available.

The effect of fluid inertia on extensional flow through constriction microchannels is characterized by the Reynolds number,

where ρ is the fluid density, V = Q/wh is the average fluid velocity in the constriction section of the channel, Q is the volumetric flow rate, w is the width of the constriction, h is the channel depth, D_h = 2wh/(w + h) is the hydraulic diameter of the constriction, is the shear-rate dependent fluid viscosity estimated at the characteristic shear rate, , in the width direction of the constriction (see in Fig. 1(b) the value of for 1 ml h⁻¹ flow rate in the 10 : 1 constriction channel). The effect of fluid elasticity is characterized by the Weissenberg number, Wi,⁴⁵

The relative importance of fluid elasticity over inertia is measured by the elasticity number, El,

which is independent of flow kinematics except for the term for shear thinning fluids. However, once the η_∞ plateaus are reached at flow rates higher than 1 ml h⁻¹ in Fig. 1(b), El remains effectively constant regardless of and can be treated as material constants thereafter. In other words, the value of Wi becomes essentially proportional to that of Re. Further, as Wi ∼ 0 in the XG solution regardless of flow rate, we will use Re to present the experimental results in the next section. The El values for the three polymer solutions at Q = 10 ml h⁻¹ in the CR = 10 channel are presented in Table 2. Note that at a fixed Q, increasing CR leads to increased Re, Wi and El because of the decreasing constriction width, w, in eqn (2)–(4).

The numerical simulation of the flow of Newtonian DI water in the three constriction microchannels was performed in COMSOL® with the Laminar Flow Module, which uses finite element method to solve the Navier–Stokes equations. The physics informed meshing option was employed to discretize the domain finely through Delauney tetrahedralization. The Dirichlet type boundary condition of no slip and no penetration were imposed on the walls. Moreover, a fully developed field was set at the inlet, and no backflow was set at the outlet. The images of predicted streamlines were extracted from the middle plane along the channel depth.

3 Results and discussion

3.1 Newtonian DI water

The contraction flow of DI water remains undisturbed, i.e., steady vortex-free, in the constriction microchannels with varying contraction ratios, CR, for flow rates up to 120 ml h⁻¹. In contrast, as viewed from the snapshot images in Fig. 2(a), the expansion flows of DI water exhibit similar feature evolution patterns among the three channels with increasing flow rates. Specifically, fluid inertia induces flow circulations in each channel above certain respective flow rates. These vortices remain symmetrical and grow along the channel length as the flow rate further increases. They start forming at the lips of the reentrant corners at the flow rate of 10 ml h⁻¹ in the CR = 5 constriction channel with Re = 37.0, whereas such onset occurs at a lower flow rate of 5 ml h⁻¹ in the CR = 10 channel with Re = 27.8. A further reduced flow rate of 3 ml h⁻¹ is observed for the formation of lip vortices in the CR = 16 channel with Re = 20.5. This tendency to form the lip vortices at lower Re with increasing CR can be explained by the transition from bimodal (with double off-center maxima) to parabolic (with a single-centered maximum) velocity profiles at the contraction, which facilitates the flow separation and vortex growth via decreasing the out-of-plane as well as spanwise momentum diffusion.³⁰ Increasing the flow rate extends the lip vortices to the salient corners of each channel and afterwards enhances the vortex size downstream with the symmetry conserved. Fig. 2(b) shows the phase diagram of flow states in the dimensionless CR−Re space. The threshold Reynolds number, Re_th, for the transition from undisturbed to vortical flows decreases with increasing CR, indicating that the inertial expansion flow of Newtonian water gets destabilized in a larger CR channel. The value of Re_th for the onset of inertial vortices in each constriction channel will be compared with those for the onset of flow instability in the three polymer solutions later in Section 3.5.

Fig. 2 Expansion flow (from top to bottom) of Newtonian water in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles, where the highlighted dimension, L_e, indicates the measurement of expansion-flow vortex length, the right-halves at 20 ml h⁻¹ are numerically simulated streamlines, and the dashed-line circles mark the lip vortices; (b) summary of the flow states in the CR−Re space, where ‘circles’ stand for undisturbed flows, ‘triangles’ for lip vortices, ‘squares’ for symmetric corner vortices, and the arrow indicates that the onset of fluid inertia induced expansion-flow vortices takes place at a smaller Re with the increase of CR; (c1) plot of the normalized vortex length, χ_e = L_e/W, against Re, where W is the width of the channel, the dashed lines are the trendlines fitted to the corresponding experimental data, and the solid lines are the numerical predictions in corresponding channels; (c2) plot of the CR modified vortex length, χ_e/CR, against Re, where the dashed line is a power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels and the solid line is from the numerical simulation.

Fig. 2(c1) compares the normalized vortex lengths, χ_e = L_e/W, as a function of Re in the three constriction channels, where L_e is the expansion flow vortex length measured directly from the images in Fig. 2(a) and W is the channel width fixed at 500 μm. While the inertial vortices are visually larger in the channel with a larger CR, the growth rate of χ_e with Re does not differ significantly among the three channels. Changing the ordinate from χ_e to χ_e/CR causes a collapse of this CR dependence as depicted in Fig. 2(c2). Such a geometric scaling law was previously shown for the abrupt contraction flow of viscoelastic fluids by Nguyen & Boger²¹ as well as Alves and co-workers,^28,32 albeit no such account in literature are known for Newtonian flows. We highlight that the vortex size collapse is experimentally robust and physically meaningful due to the simplicity of the modified parameter, which conveys a linear CR dependence of the vortex size within the range of planar constriction channel dimensions considered here. To verify the inertial feature of such expansion flow vortices, we conducted 3D numerical simulations of the Newtonian flow field in each constriction channel for a range of flow rates. As seen from Fig. 2(a), the predicted streamlines for the flow rate of 20 ml h⁻¹ appear to closely match the experimental pathlines of tracer particles in each constriction channel. The predicted vortex lengths, χ_e, vs. Re in each channel also show a good agreement with the experimental data in Fig. 2(c1). Moreover, as demonstrated in Fig. 2(c2), the predicted CR modified vortex lengths, χ_e/CR, from all three channels also collapse to a power-law line that nearly overlaps with the trendline fitted to the experimental data. The discrepancy may stem from the differences in the experimental and the numerical conditions such as the channel geometry, dimensions, corner rounding, surface roughness, and deformation of PDMS under higher strains.⁴⁶

3.2 Shear thinning XG solution

Fig. 3(a) shows the snapshots of the contraction flow of XG solution in the three constriction microchannels. The lowest tested flow rate is 0.1 ml h⁻¹ with Re ≤ 0.11, at which 2D symmetric vortices spanning from the constriction lips to the salient corners are observed to have formed in all three channels. These fluid shear thinning-induced vortices grow in length with increasing flow rates and start overlapping with apparent out-of-plane circulations at 20 ml h⁻¹ in the CR = 5 channel with Re = 30.9 because of the coupled impact of fluid inertia and shear thinning. They remain quasistatic in time with small fluctuations in the boundaries compared to their time-averaged sizes. Further increasing the flow rate results in an enlargement of these quasistatic vortices, which, however, remain separated without any direct interactions. The transition to the quasistatic state takes place at 10 ml h⁻¹ in the CR = 10 channel with Re = 25.5 compared to 3 ml h⁻¹ in the CR = 16 channel with Re = 9.11. These flow states are summarized in the CR−Re space in Fig. 3(b), where it is evident that the flow transition to instability is incurred monotonically at lower Re as the channel CR is increased. Such trend perhaps could be explained via the evolution of velocity profile with CR and the associated off-plane and spanwise momentum transfer in these shear-thinning flows as Oliveira et al.³⁰ did for Newtonian flows. Fig. 3(c) plots the normalized vortex length, χ_c = L_c/w, against Re, where L_c is the contraction-flow vortex length measured directly from the images in Fig. 3(a). The vortex sizes coincide along a single line regardless of CR, which, as highlighted in Fig. 3(c), follows a power dependence on Re. A similar geometric independence of the contraction-flow vortex size from the constriction length was observed in the XG solution in our previous study.⁴⁰ These observations may imply the competitive effects of fluid shear thinning (enhancement) and inertia (suppression) in the contraction part of the channel.

Fig. 3 Contraction flow (from top to bottom) of shear thinning XG solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles, where the highlighted dimension, L_c, indicates the measurement of contraction-flow vortex length; (b) summary of the flow states in the CR−Re space, where ‘squares’ stand for 2D corner vortices, ‘diamonds’ for 3D overlapping corner vortices, and the arrow indicates that the onset of overlapping vortices takes place at a smaller Re with the increase of CR; (c) plot of the normalized vortex length, χ_c = L_c/W, against Re, where the dashed line is the power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels.

Fig. 4(a) shows the snapshots of XG expansion flow in the three constriction channels. For the CR = 10 channel, the streamlines near the constriction lips start deviating from smooth paths of undisturbed state to a weakly unsteady bending state at 10 ml h⁻¹ (Re = 25.5). Spike-like secondary flow zones appear on the constriction lips at 20 ml h⁻¹, where two mildly fluctuating vortices form simultaneously at the salient corners. These pairs of spikes and salient corner vortices maintain symmetry about the channel centerline. Increasing the flow rate enlarges these vortices causing them to merge at 30 ml h⁻¹ (Re = 84.1) on their respective sides, which maintain the spatial pseudo-symmetry (i.e., the quasistatic vortex on one side of the expansion remains similar in time-averaged size to that on the other side) and grow in length with increasing flow rates. Decreasing CR to 5 demonstrates an overall similar flow feature pattern to that in the 10 : 1 channel, whereas the transition to streamline bending occurs at a higher flow rate of 30 ml h⁻¹ (Re = 49.2). On the other hand, increasing CR to 16 incurs a similar sequential development of states as the other two channels, albeit earlier transitions can be noticed. The evolution of flow states in each of the constriction channels is displayed in the CR−Re space in Fig. 4(b). Consistent with the contraction flow in Fig. 3(b), the threshold Re for the transition from smooth to bending streamlines in the expansion flow also gets diminished with increasing CR. Fig. 4(c1) shows the normalized expansion flow vortex length, χ_e, which follows a power-law dependence on Re in each channel and gets larger with the increase of CR for the same Re. Fig. 4(c2) shows the CR modified expansion vortex length, χ_e/CR, vs. Re, where the experimental data from all three channels collapse into one power trendline with an index of almost 2.

Fig. 4 Expansion flow (from top to bottom) of shear thinning XG solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles illustrate the streamlines; (b) summary of the flow states in the CR−Re space, where ‘circles’ stand for undisturbed flows, ‘rectangles’ for bending streamlines, ‘triangles’ for lip vortices, ‘squares’ for salient corner vortices, and the arrow indicates that the onset of flow instability takes place at a smaller Re with the increase of CR; (c1) plot of the normalized vortex length, χ_e, against Re, where the dashed lines are the trendlines fitted to the corresponding experimental data; (c2) plot of the CR modified vortex length, χ_e/CR, against Re, where the dashed line is a power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels.

3.3 Viscoelastic PEO solution

Fig. 5(a) shows the snapshots of the contraction flow of PEO solution, where the feature evolution with flow rate overall follows the same pattern in all three channels. However, the feature persistence within flow rate intervals varies based on the channel CR. The flow in CR = 10 channel transitions from undisturbed to bending streamlines at 8 ml h⁻¹ (Re = 20.4). A single vortex develops near the salient corner of the contraction at around 20 ml h⁻¹ and grows in circulation size and intensity till the highest tested flow rate of 50 ml h⁻¹ (Re = 127). The CR = 5 channel sees a flow transition to bending streamlines at 10 ml h⁻¹ (Re = 17.0) and single vortex at 35 ml h⁻¹. However, the latter state persists only till 40 ml h⁻¹ (Re = 68.0), when a pair of similar sized vortices appear on each side of the contraction. These double vortices are loosely stable in time as they mildly pulsate in the longitudinal direction while maintaining a pseudo symmetry. They grow in length without altering the state for flow rates up to 120 ml h⁻¹ (Re = 204). In the CR = 16 channel, bending streamlines show up at 3 ml h⁻¹ (Re = 9.41). This state directly transitions to the pseudo symmetric double vortex state at 8 ml h⁻¹ (Re = 25.1) bypassing the single vortex state. The vortex size enlargement continues till 30 ml h⁻¹ (Re = 94.1) when the symmetry significantly breaks as one of the vortices becomes much longer and wider than the other. These asymmetric vortices demonstrate bistable properties, suddenly exchanging sides at random intervals.

Fig. 5 Contraction flow (from top to bottom) of viscoelastic PEO solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles; (b) summary of the flow states in the CR−Re space, where ‘circles’ stand for undisturbed flow, ‘rectangles’ for bending streamlines, ‘diamonds’ for pseudo-symmetric double vortices, ‘pluses’ for single vortex state, ‘crosses’ for asymmetric double vortices, and the two arrows indicate that the threshold Re for the onset of flow instability exhibits a nonmonotonic dependence on CR; (c) plot of the normalized vortex length, χ_c, against Re, where the dashed lines are the logarithmic trendlines fitted to the corresponding experimental data in each channel.

Fig. 5(b) summarizes the flow states, where the threshold Re for the onset of bending streamlines first increases and then decreases with increasing CR, exhibiting a surprising nonmonotonic trend. The comparison of the vortex size is shown in Fig. 5(c), where a nonmonotonic dependence of χ_c on the channel CR is also observed with the vortices in the CR = 10 channel being the smallest. Such nonmonotonic behaviors are absent in our earlier experiments of the constriction length⁴⁰ and depth⁴¹ effects on the same PEO solution flow. They were not observed either in previous studies that were mostly concerned with creeping flows of viscoelastic fluids.^5,12 We note that significant inertia is involved in the non-vortical to vortical state transition in the viscoelastic PEO flow, whereas such transition takes place under creeping conditions in the shear-thinning XG flow that observes a monotonic trend. Hence, we hypothesize that the nonmonotonic trend might be a signature nonlinear effect in elasto-inertial flows and could be better explained with insights on the change of velocity profile and momentum transfer due to CR in such flows at the contraction. Such trend may also imply that the favorability of polymer conformational transition follows a nonmonotonic dependence on CR as polymers pass through the constriction region.

Fig. 6(a) shows the expansion flow snapshots of PEO solution in the three constriction channels. In the CR = 5 channel, the flow transitions from undisturbed to an unsteady state of streamlines bending convergently towards the center at around 15 ml h⁻¹ (Re = 25.5). Then at 20 ml h⁻¹ (Re = 34.0), a pair of small secondary flow circulations appear at the two salient corners of the expansion. They sluggishly grow in length and width with increasing flow rates to 120 ml h⁻¹ (Re = 204). In the CR = 10 channel, the bending of streamlines starts earlier at 8 ml h⁻¹ (Re = 20.4), so does the salient corner vortices. In the CR = 16 channel, the transitions to streamline bending and vortex onset occur at 5 ml h⁻¹ (Re = 15.7) and 10 ml h⁻¹ (Re = 31.4), respectively. The vortices again initiate from the salient corners but grow significantly with increasing flow rates. They even reach the constriction lips from the salient corners in width at flow rates above 20 ml h⁻¹. This facilitation of vortices in the CR = 16 channel in contrast to their suppression in the CR = 5 channel could be occurring due to possibly enhanced polymer degradation in the slit-like former channel. Recent works by Yin et al.^11,19,39 have revealed that PEO with similar molecular properties can undergo severe scissions and thus experience reduced elasticity while flowing through such channels, which bolsters the possibilities of similar occurrences in our flows. Fig. 6(b) summarizes the flow states in the CR−Re space, where the threshold Re for the onset of bending streamlines decreases with increasing CR like in the expansion flow of all other fluids tested in this work. Fig. 6(c1) shows the normalized vortex length χ_e against Re, which is noticeably greater in the larger CR channel while its growth rate does not differ significantly. Fig. 6(c2) plots the CR modified vortex length, χ_e/CR, against Re, where the data points from all three channels again collapse into a power trendline with an index of approximately 1, smaller than 1.4 in water and 2 in the XG solution. This declined vortex growth rate in the PEO solution may be attributed to the viscoelastic suppression of inertial flow separation that were demonstrated in our previous studies.^40,41,47,48

Fig. 6 Expansion flow (from top to bottom) of viscoelastic PEO solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles illustrate the streamlines; (b) summary of the flow states in the CR−Re space, where ‘circles’ stand for undisturbed flows, ‘rectangles’ for bending streamlines, ‘triangles’ for small vortices at the salient corners, ‘squares’ for full salient corner vortices, and the arrow indicates that the onset of streamline bending takes place at a smaller Re with the increase of CR; (c1) plot of the normalized vortex length, χ_e, against Re, where the dashed lines are the trendlines fitted to the corresponding experimental data; (c2) plot of the CR modified vortex length, χ_e/CR, against Re, where the dashed line is a power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels.

3.4 Shear thinning and viscoelastic PAA solution

Fig. 7(a) shows the contraction flow snapshots of the PAA solution. In the CR = 10 channel, large steady symmetric vortices form at the salient corners of the contraction at 0.1 ml h⁻¹ (Re = 0.170) because of the strong fluid shear thinning effect. They grow in length with increasing flow rates till 2 ml h⁻¹ (Re = 5.00) when the vortices become asymmetric and begin showing aperiodic spanwise oscillations. Then, the flow becomes unsteady at 10 ml h⁻¹ (Re = 26.0), which gets enhanced visually in intensity for the rest of the flow rates tested. This feature pattern is identical to our earlier observations on the PAA flow through similar constriction channels.^40,41,47 Increasing CR to 16 also maintains the same pattern except that the flow transitions occur at lower corresponding flow rates. Surprisingly, the flow pattern is drastically different in the CR = 5 channel as the vortices transition from steady symmetric to weakly fluctuating yet spatially symmetric at 10 ml h⁻¹ (Re = 16.8). This state visually resembles the quasi-static vortex pairs observed in the contraction flow of XG solution at higher inertias in Fig. 3(a). The state of asymmetric fluctuating vortices does not incur till the highest tested flow rate of 120 ml h⁻¹ in this channel with Re = 210. A summary of the flow states is presented in Fig. 7(b), where the threshold Re for the onset of asymmetric fluctuating vortices decreases with increasing CR. Fig. 7(c1) plots the normalized vortex length, χ_c, which increases in a larger CR channel at the same Re. This monotonic dependence is different from the independence and the nonmonotonic dependence of χ_c on CR in the contraction flows of XG and PEO solutions, respectively, because of probably the simultaneous shear thinning and viscoelastic features of the PAA solution. Fig. 7(c2) plots the CR modified vortex length, , against Re, where the data points from all three channels collapse into one power trendline with an index of about 0.2.

Fig. 7 Contraction flow (from top to bottom) of shear thinning and viscoelastic PAA solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles; (b) summary of the flow states in the CR−Re space, where ‘squares’ stand for steady symmetric corner vortices, ‘diamonds’ for pseudo-symmetric double vortices, ‘crosses’ for asymmetric fluctuating vortices, ‘asterisks’ for unsteady flow, and the arrow indicates that the threshold Re for the onset of asymmetric fluctuating vortices decreases with the increase of CR; (c1) plot of the normalized vortex length, χ_c, against Re, where the dashed lines are the trendlines fitted to the corresponding experimental data; (c2) plot of the CR modified vortex length, , against Re, where the dashed line is a power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels.

Fig. 8(a) shows the expansion flow snapshots of the PAA solution. Overall, the flow states in each constriction channel follow a general trend of transitioning from undisturbed to streamline bend and then the formation of a pair of quasi-static lip spikes and salient corner vortices on both sides of the expansion. Following that are the spike-vortex merging and pseudo-symmetric quasi-static vortex enhancement regimes with flow rate increments. The differences among the three channels are that the transition to every flow state takes place at a smaller Re in a larger CR channel. For example, the spike formation initiates at 40 ml h⁻¹ in the 5 : 1 channel with Re = 69.3, at 15 ml h⁻¹ in the 10 : 1 channel with Re = 39.2, and at 8 ml h⁻¹ in the 16 : 1 channel with Re = 32.3. Fig. 8(b) shows the flow state summary in the CR−Re space, which clearly demonstrates the flow destabilization in a larger CR channel. Fig. 8(c1) plots the normalized vortex length, χ_e, which exhibits a positive power-law dependence on the channel CR and becomes the largest in the 16 : 1 channel for the same values of Re. Fig. 8(c2) replots this vortex length in the χ_e/CR²−Re space, where the data points from all three channels collapse into a power trendline. This quadratic dependence of χ_e on CR in the PAA expansion flow is stronger than the linear relationship in the expansion flow of all other tested fluids. More studies, especially numerical,⁴⁹ will be needed to understand our experimentally observed variations of contraction and expansion flow patterns among fluids with distinct rheological properties.

Fig. 8 Expansion flow (from top to bottom) of shear thinning and viscoelastic PAA solution in constriction microchannels with three contraction ratios, CR, at varying flow rates: (a) snapshot images of tracer particles; (b) summary of the flow states in the CR−Re space, where ‘circles’ stand for undisturbed flows, ‘rectangles’ for bending streamlines, ‘triangles’ for coexisting lip and salient corner vortices, ‘squares’ for salient corner vortices, and the arrow indicates that the onset of streamline bending takes place at a smaller Re with the increase of CR; (c1) plot of the normalized vortex length, χ_e, against Re, where the dashed lines are the trendlines fitted to the corresponding experimental data; (c2) plot of the CR modified vortex length, χ_e/CR², against Re, where the line is a power trendline [equation (x and y represent the variables on the horizontal and vertical axes, respectively) and R² value highlighted on chart were obtained from Microsoft Excel®] fitted to the experimental data from all three channels.

3.5 Comparison among the tested fluids

Fig. 9(a) shows the variations of the threshold Reynolds number, Re_th, against the channel CR for the onset of contraction flow instability in XG (overlapping vortices), PEO (bending streamlines), and PAA (fluctuating vortices) solutions, respectively. Increasing CR leads to a reduced Re_th in the shear thinning XG and PAA solutions. In contrast, the value of Re_th in the viscoelastic PEO solution first increases and then decreases with the increase of CR from 5 to 10 and then 16. Such a nonmonotonic dependence on CR also applies to the normalized contraction flow vortex length, χ_c, in the PEO solution, which, as shown in Fig. 5(c), has the smallest value in the CR = 10 channel for the same Re. It is, however, absent from the two shear thinning fluids, where χ_c is independent of CR in the XG solution while having a weak 1/3-order dependence on CR in the PAA solution. The latter index value can be viewed as the consequence of the coupled influence of shear thinning and elasticity in the PAA solution. These phenomena together may indicate the opposite and stronger impact of shear thinning than elasticity on the contraction flow, consistent with the observations in our recent studies on the effects of constriction length and depth, respectively, in microfluidic extensional flow.^40,41

Fig. 9 Comparison of the threshold Reynolds number, Re_th, for the onset of flow instability in the tested fluids through constriction microchannels with varying CR: (a) Re_th in the contraction flow; (b1) Re_th and (b2) in the expansion flow. All the lines are used to guide the eyes only.

Fig. 9(b1) demonstrates that increasing CR continuously lowers Re_th for the onset of first state change in the expansion flow of tested fluids. Specifically, the first disturbed state following the undisturbed flow in water is the appearance of fluid inertia-induced lip vortices, whereas it is the bending of streamlines in the three polymer solutions. The values of Re_th in water are roughly equal to those in the XG solution for all three channels. They are, however, apparently higher than those in the PEO solution and reach the lowest in the PAA solution in each constriction channel. These phenomena indicate the dominant and suppression effect of fluid elasticity over shear thinning on the expansion flow instability, again consistent with our recent studies of the constriction length and depth effects.^40,41Fig. 9(b2) plots against CR, which shows an interesting self-similarity in the trends for all tested fluids. In other words, the first disturbed expansion flow state initiates at a nearly constant regardless of the channel CR in every tested fluid. A qualitatively similar trend was noted in previous studies from Alves and co-workers,^28,32 where the modification of the Deborah number, De, by CR to De/CR leads to the onset similarities in viscoelastic contraction flows. This self-similarity phenomenon is not observed from the plots of normalized expansion flow vortex length, χ_e, which scales linearly with CR in all fluids but the PAA solution (quadratically) for the same Re. Moreover, the CR modified χ_e follows a power-law dependence on Re in each fluid but with distinct index values.

It is imperative to note that the disparate trends and features in the three fluids may not be fully described with conventional rheometric characterization alone. The polymer chain dynamics and interactions in microchannel confinements can often incur stronger non-Newtonian responses that may not be discernible in macroscopic rheometers. For example, blood plasma was elusively characterized for decades as Newtonian with conventional rheometers. The notion was overturned with experimental evidence for viscoelasticity in contraction-expansion microchannels through extra pressure drop measurements, which was later numerically validated for such flows.^50,51 In light of such implications, our current study unveils the characters of these polymers of practical importances in micro- and macro- environments informing crucial scale dependent properties of their flows. In addition, the migration of polymer chains giving rise to nonhomogeneous distribution and local rheology may have impacted our observed flow patterns consequentially. For instance, a Boger fluid like PEO may encounter a strong depletion of polymers near the walls for large constriction ratios.⁵² Large molecules with high rigidity, such as DNA, tend to move faster towards the center in microchannels with higher intrachain interactions.^53,54 However, the existing works have been conducted mostly under creeping conditions, which can have highly disparate trends in contrast to considerably inertial flows. Insights on how CR impacts the migration of polymers with different rigidity, size and polydispersity, such as those used in this work, under significant inertia are yet to be unearthed. Hence, it is not feasible to draw any systematic or robust conclusions based on polymer migrations for our observed flow feature developments. Nevertheless, we note that significant contributions from the migration phenomena in polymer solutions may be expected in our flow fields.

4 Conclusions

We have experimentally studied the effects of changing contraction ratio, CR, on the flow of rheologically distinct polymer solutions compared to Newtonian water in planar constriction microchannels. Flow visualizations reveal that increasing the channel CR within the limit of our tests enhances the expansion flow vortices monotonically in each fluid. Moreover, the vortex size scales linearly with CR in all fluids but the PAA solution, where a quadratic dependence on CR is observed because of perhaps its simultaneously shear thinning and elastic features. The threshold Reynolds number, Re_th, for the initiation of expansion flow disturbances acts as a monotonic negative function of CR in each fluid. This dependence, however, vanishes if the disturbance onsets are expressed in terms of , which then becomes solely a function of fluid rheology. On the other hand, the contraction flow vortex size has a zero or 1/3 order dependence on CR in strongly shear thinning XG and PAA flows. In contrast, the contraction flow vortices in the purely viscoelastic PEO solution scale nonmonotonically with CR, first diminishing and then rising. The values of Re_th for the onset of contraction flow instability in the PEO solution also follow a nonmonotonic correlation with CR, whereas those in the shear thinning XG and PAA solutions both get smaller with increasing CR. We hope that the results presented here will stimulate further experimental and theoretical (including numerical) studies on similar topics, such as transient nonlinear rheological characterization of test fluids and tube theory-based modeling of contraction-expansion flows,⁵⁵ for improved understanding of polymer solution flow transitions and scaling laws.

Author contributions

M. R. performed the experiment and conducted the analysis; M. M., D. H., Y. S. and X. X. designed the project; X. X. supervised the project; M. R. wrote the manuscript; M. M., D. H. and X. X. edited the manuscript; all authors commented on the manuscript.

Conflicts of interest

None.

Data availability

The data supporting this article are within the article.

Acknowledgements

This work was supported in part by Clemson University through a SEED grant (X. X.) and University 111 Project of China under grant number B08046 (Y. S.).

References

1. C. Yuan, H.-N. Zhang, Y.-K. Li, X.-B. Li, J. Wu and F.-C. Li, Proc. Inst. Mech. Eng., Part C, 2020, 234, 4390–4414,  DOI:10.1177/0954406220922863.
2. A. C. Barbati, J. Desroches, A. Robisson and G. H. McKinley, Annu. Rev. Chem. Biomol. Eng., 2016, 7, 415–453.
3. X. Lu, C. Liu, G. Hu and X. Xuan, J. Colloid Interface Sci., 2017, 500, 182–201.
4. A. Anbari, H. T. Chien, S. S. Datta, W. Deng, D. A. Weitz and J. Fan, Small, 2018, 14, 1703575.
5. S. S. Datta, A. M. Ardekani, P. E. Arratia, A. N. Beris, I. Bischofberger, G. H. McKinley, J. G. Eggers, J. E. López-Aguilar, S. M. Fielding, A. Frishman, M. D. Graham, J. S. Guasto, S. J. Haward, A. Q. Shen, S. Hormozi, A. Morozov, R. J. Poole, V. Shankar, E. S. G. Shaqfeh, H. Stark, V. Steinberg, G. Subramanian and H. A. Stone, Perspectives on viscoelastic flow instabilities and elastic turbulence, Phys. Rev. Fluids, 2022, 7, 080701,  DOI:10.1103/PhysRevFluids.7.080701.
6. B. Wu, S. Liu, D. Jiang and W. Tang, Electrophoresis, 2024, 45, 1233–1242.
7. L. E. Rodd, J. J. Cooper-White, D. V. Boger and G. H. McKinley, J. Non-Newtonian Fluid Mech., 2007, 143, 170–191.
8. F. J. Galindo-Rosales, L. Campo-Deaño, P. C. Sousa, V. M. Ribeiro, M. S. N. Oliveira, M. A. Alves and F. T. Pinho, Exp. Therm. Fluid Sci., 2014, 59, 128–139.
9. C. A. Browne, A. Shih and S. S. Datta, Small, 2020, 16, 1903944.
10. F. Hillebrand, S. Varchanis, C. C. Hopkins, S. J. Haward and A. Q. Shen, Soft Matter, 2024, 20, 7133–7146.
11. G. Yin, Y. Nakamura, H. Suzuki, F. Lequeux and R. Hidema, Phys. Fluids, 2024, 36, 043114.
12. C. Sasmal, J. Non-Newtonian Fluid Mech., 2025, 337, 105393.
13. G. H. McKinley, W. P. Raiford, R. A. Brown and R. C. Armstrong, J. Fluid Mech., 1991, 223, 411–456.
14. M. A. Alves and R. J. Poole, J. Non-Newtonian Fluid Mech., 2007, 144, 140–148.
15. L. E. Rodd, D. Lee, K. H. Ahn and J. J. Cooper-White, J. Non-Newtonian Fluid Mech., 2010, 165, 1189–1203.
16. C. A. Browne, A. Shih and S. S. Datta, J. Fluid Mech., 2020, 890, A2.
17. D. W. Carlson, A. Q. Shen and S. J. Haward, J. Fluid Mech., 2021, 923, R6,  DOI:10.1017/jfm.2021.620.
18. C. C. Hopkins, S. J. Haward and A. Q. Shen, Soft Matter, 2022, 18, 4868–4880.
19. G. Yin, Y. Nakamura, H. Suzuki, F. Lequeux and R. Hidema, Phys. Fluids, 2024, 36, 124139.
20. E. Y. Chen and S. S. Datta, J. Rheol., 2025, 69, 235–249.
21. H. Nguyen and D. V. Boger, J. Non-Newtonian Fluid Mech., 1979, 5, 353–368.
22. R. E. Evans and K. Walters, J. Non-Newtonian Fluid Mech., 1986, 20, 11–29.
23. R. E. Evans and K. Walters, J. Non-Newtonian Fluid Mech., 1989, 32, 95–105.
24. K. Chiba, T. Sakatani and K. Nakamura, J. Non-Newtonian Fluid Mech., 1990, 36, 193–203.
25. G. H. McKinley, W. P. Raiford, R. A. Brown and R. C. Armstrong, J. Fluid Mech., 1991, 223, 411–456.
26. J. M. Marchal and M. J. Crochet, J. Non-Newtonian Fluid Mech., 1987, 26, 77–114.
27. B. Purnode and M. J. Crochet, J. Non-Newtonian Fluid Mech., 1996, 65, 269–289.
28. M. A. Alves, P. J. Oliveira and F. T. Pinho, J. Non-Newtonian Fluid Mech., 2004, 122, 117–130.
29. M. S. N. Oliveira, M. A. Alves, H. G. McKinley and F. T. Pinho, 13th International Symposium on Applications of Laser Techniques to Fluid Mechanics, Lisbon, Portugal, 26-29 June, 2006.
30. M. S. N. Oliveira, M. A. Alves, F. T. Pinho and G. H. McKinley, Exp. Fluids, 2007, 43, 437–451.
31. P. C. Sousa, I. S. Pinho, F. T. Pinho, M. S. N. Oliveira and M. A. Alves, Comput. Methods Appl. Sci., 2011, 19, 265–279.
32. M. S. N. Oliveira, P. J. Oliveira, F. T. Pinho and M. A. Alves, J. Non-Newtonian Fluid Mech., 2007, 147, 92–108.
33. P. C. Sousa, P. M. Coelho, M. S. N. Oliveira and M. A. Alves, J. Non-Newtonian Fluid Mech., 2009, 160, 122–139.
34. P. C. Sousa, P. M. Coelho, M. S. N. Oliveira and M. A. Alves, Chem. Eng. Sci., 2011, 66, 998–1009.
35. A. Lanzaro and X. F. Yuan, J. Non-Newtonian Fluid Mech., 2011, 166, 1064–1075.
36. M. Pérez-Camacho, J. E. López-Aguilar, F. Calderas, O. Manero and M. F. Webster, J. Non-Newtonian Fluid Mech., 2015, 222, 260–271.
37. J. E. López-Aguilar, M. F. Webster, H. R. Tamaddon-Jahromi, M. Pérez-Camacho and O. Manero, J. Non-Newtonian Fluid Mech., 2016, 237, 39–53.
38. R. Hidema, T. Oka, Y. Komoda and H. Suzuki, Phys. Fluids, 2019, 31, 072005.
39. G. Yin, Y. Nakamura, H. Suzuki, F. Lequeux and R. Hidema, Phys. Fluids, 2024, 36, 124144.
40. M. K. Raihan, S. Wu, Y. Song and X. Xuan, Soft Matter, 2021, 17, 9198–9209.
41. M. K. Raihan, S. Wu, H. Dort, M. Baghdady, Y. Song and X. Xuan, Soft Matter, 2022, 18, 7427–7440.
42. A. Lindner, D. Bonn and J. Meunier, Phys. Fluids, 2000, 12, 256–261.
43. L. E. Rodd, T. P. Scott, D. V. Boger, J. J. Cooper-White and G. H. McKinley, J. Non-Newtonian Fluid Mech., 2005, 129, 1–22.
44. R. J. Poole and M. P. Escudier, J. Non-Newtonian Fluid Mech., 2004, 117, 25–46.
45. R. B. Bird, R. C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Wiley-Interscience, Hoboken, NJ, USA, 1987, vol. 1.
46. D. Li, L. Song, C. Zhang, L. Yu and X. Xuan, Phys. Fluids, 2021, 33, 012002.
47. P. P. Jagdale, D. Li, X. Shao, J. B. Bostwick and X. Xuan, Micromachines, 2020, 11, 278.
48. S. Wu, M. K. Raihan, L. Song, X. Shao, J. B. Bostwick, L. Yu, X. Pan and X. Xuan, J. Non-Newtonian Fluid Mech., 2021, 290, 104508.
49. M. A. Alves, P. J. Oliveira and F. T. Pinho, Annu. Rev. Fluid Mech., 2021, 53, 509–541.
50. M. Brust, C. Schaefer, R. Doerr, L. Pan, M. Garcia, P. E. Arratia and C. Wagner, Phys. Rev. Lett., 2013, 110, 078305.
51. S. Varchanis, Y. Dimakopoulos, C. Wagner and J. Tsamopoulos, Soft Matter, 2018, 14, 4238–4251.
52. S. Tsouka, Y. Dimakopoulos, V. Mavrantzas and J. Tsamopoulos, J. Rheol., 2014, 58, 911–947.
53. K. A. Dill and B. H. Zimm, Nucleic Acids Res., 1979, 7, 735–749.
54. H. Ma and M. D. Graham, Phys. Fluids, 2005, 17, 083103.
55. S. Varchanis, Y. Dimakopoulos and J. Tsamopoulos, J. Rheol., 2018, 62, 25–47.

---