Tristability in viscoelastic flow past side-by-side microcylinders

Since the advent of microfluidics in the early 2000s [1, 2], geometries with length-scales $\ell\sim O(100~{}\mu$m) have become a vital tool in experimental fluid dynamics. At the microscale, viscoelastic fluids (with properties between viscous liquids and elastic solids) can flow with negligible inertia (Reynolds number $\textnormal{Re}\sim\ell\ll 1$), but high elasticity (Weissenberg number $\textnormal{Wi}\sim\ell^{-1}\gg 1$) [2]. In such flows, elasticity becomes the dominant source of nonlinearity, leading to instabilities [3, 4, 5, 6, 7, 8, 9], and time-dependency that impact widespread processes ranging from jet fragmentation [10, 11] to hemodynamics [12, 13] and porous media flows [14, 15, 16, 17, 18, 19, 20]. In particular, the path selection and switching phenomena in viscoelastic porous media flow is considered of fundamental importance in processes such as enhanced oil recovery, groundwater remediation, filtration, and drug delivery  [14, 15, 16, 17, 18, 19].

Porous media are frequently modeled by ordered and disordered arrays of microfluidic circular cylinders [15, 16, 17, 18, 19]. Flow past a single circular cylinder in a channel is an archetypal problem in fluid dynamics, and a ‘benchmark’ for studying viscoelastic flows. The stagnation point downstream of a cylinder is a location where streamline curvature combines with strong velocity gradients; conditions that render viscoelastic base flows prone to instability and downstream fluctuations [21, 22, 23, 24]. For fluids with a shear-rate-dependent viscosity (i.e., shear-thinning), the perturbation to the base flow can lead to a steady symmetry-breaking flow bifurcation where the viscoelastic fluid selects a preferred path around one, or other, side of the cylinder [5, 6, 7, 25]. This behavior has clear relevance to understanding transport through porous arrays, but the interaction with neighboring array elements is lacking. Building ‘bottom-up’ complexity towards more realistic model systems, it is natural to consider two cylindrical objects either aligned in the flow direction, or positioned side-by-side in a channel. Viscoelastic flow past two (or more) objects aligned on the flow axis is a well studied problem (e.g., Refs. 23, 26, 27, 8). However, although equally important, the case of two objects positioned transverse to the flow has received scant attention, with only one numerical study conducted at high Reynolds number [28]. To date, creeping viscoelastic flow past side-by-side cylinders has not been studied.

In this Letter, we present microfluidic experiments of a viscoelastic shear-thinning fluid flowing past two microcylinders transverse to the primary flow direction (Fig. 1) and show that the resulting non-linear flow behavior at high Wi is significantly influenced by the spacing of the cylinders. We show that, due to a combination of supercritical bifurcations that occur as Wi is varied, multiple stable flow states are possible in a given geometry. This is the first study of low-Re viscoelastic flow in such a geometry and serves as a fundamental contribution towards understanding deterministic path-selection in porous media flow.

The model viscoelastic fluid is a well-studied aqueous wormlike micellar (WLM) solution consisting of 100 mM cetylpyridinium chloride (CPyCl) and 60 mM sodium salicylate (NaSal) [29, 30]. At $24^{\circ}$C (ambient laboratory temperature), the entangled WLM solution has a zero shear viscosity $\eta_{0}=47$ Pa s, exhibits a stress-plateau (shear-banding region [31]), and in small-amplitude oscillation is well-described by a single-mode Maxwell model with relaxation time $\lambda=1.7$ s (Fig. S1 [32]).

Microfluidic channels (Fig. 1) were fabricated in fused silica by selective laser-induced etching [33]. The eleven channels used all have a rectangular cross-section with width $W=400~{}\mu$m transverse to the flow ($y$-direction) and height $H=2000~{}\mu$m in the neutral ($z$) direction. Each channel contains two cylinders of radius $R=20~{}\mu$m equally spaced either side of the primary flow ($x$) axis. The intercylinder separation $L_{1}$ is varied between channels in the range $107<L_{1}<147~{}\mu$m. The spacing between the cylinders and the channel sidewalls is $L_{2}=(W-L_{1}-4R)/2$, and we define a dimensionless gap ratio $G=L_{1}/(L_{1}+L_{2})$. This parameter in principle spans $0<G<1$, where $G=0$ implies the two cylinders are touching at the channel centerline, while $G=1$ implies the cylinders are touching opposite channel walls. The channels used span a range $0.50\leq G\leq 0.62$, which encompasses the full range of flow behavior.

Flow is driven by syringe pumps (Cetoni GmbH) programmed to impose quasistatic variations in the volumetric flow rate $Q$, hence average flow velocity ${U=Q/WH}$, and Weissenberg number Wi $=\lambda U/R$. Quantitative spatially-resolved flow fields are obtained using micro-particle image velocimetry ($\mu$-PIV, TSI Inc., MN [34, 35]). At each imposed Wi, the motion of a low concentration of fluorescent seeding particles ($2~{}\mu$m diameter) is captured at the channel half-height ($z=0$ plane) using an inverted microscope (Nikon Ti) with a $5\times$, NA = 0.15 numerical aperture objective lens and a high speed camera (Phantom Miro) working in frame-straddling mode at 25 Hz. Cross-correlation between images yields velocity vectors ${\bf{u}}=(u,v)$. Since the flows examined are all time invariant, data are ensemble averaged over a 6 s sampling window. Due to shear-localization at the channel walls, the flow profile is essentially plug-like over most of the channel cross-section [6]. Therefore, the shear-rate near the cylinders is small and we define $\textnormal{Re}=\rho UR/\eta_{0}$, where $\rho=1000$ kg m^-3 is the fluid density. In all experiments, $\textnormal{Re}\lesssim 10^{-4}$.

Flow fields representative of those observed as Wi is varied are shown in Fig. 2 using two channels with contrasting $G$. Fig. 2a-c and Fig. 2d-f illustrate the behavior for ‘small’ and ‘large’ $G$, respectively. Irrespective of $G$, for low $\textnormal{Wi}<\textnormal{Wi}_{1}\approx 15$ (Figs. 1a,d), elastic and inertial forces are small and the flow is dominated by the viscous force. Flow is approximately symmetric about $x=0$ and $y=0$, and fluid passes through all three available gaps. For ‘small’ $G=0.500$, as Wi exceeds $\textnormal{Wi}_{1}$ (Fig. 2b), elasticity dominates and the system undergoes a first transition from the low-Wi symmetric state to a diverging ‘D’ state where the fluid avoids the gap between the cylinders and flows symmetrically around their sides. The velocity field is qualitatively similar to that for viscoelastic flow around a single obstacle [27, 6]. For a Newtonian fluid, two objects appear as one when the ratio of separation to radius, $L_{1}/R<0.2$ [36], whereas here the ratio is much greater at $L_{1}/R>5$. Further increasing Wi, the system undergoes a second transition at $\textnormal{Wi}_{2}\approx 50$ to an asymmetric-diverging ‘AD’ state in which the fluid selects a single preferred path either above ($y>0$, Fig. 2c₁), or below ($y<0$, Fig. 2c₂) the pair of cylinders. This randomly chosen bias is also similar to that observed for viscoelastic shear-thinning fluids flowing around a single cylinder [6, 7, 25].

For ‘large’ $G=0.603$, the first transition at $\textnormal{Wi}>\textnormal{Wi}_{1}$ results in a converging ‘C’ flow state where the fluid flows preferentially between the cylinders, avoiding the gaps at their sides (Fig. 2e). In contrast to the small-$G$ case, as Wi increases there is no second transition at $\textnormal{Wi}_{2}$ and the C state is maintained until the flow eventually becomes time-dependent at $\textnormal{Wi}\gg\textnormal{Wi}_{1}$. The nature of the time-dependence is qualitatively similar to that previously reported for a single cylinder [6] and will not be discussed further here. We note that results for small and large $G$ in Fig. 2 using a shear-thinning, but non-shear-banding viscoelastic polymer solution, show analogous flow behavior see Figs. S2 and S3 [32].

We quantify the critical flow behavior using two dimensionless flow asymmetry parameters $I^{\prime}$ and $I^{\prime\prime}$:

Here, $\bar{u}_{+}$, $\bar{u}_{-}$, and $\bar{u}_{0}$ are the average values of $u$ in the upper, lower, and intercylinder gaps, respectively (see Fig. 1). $I^{\prime}$ serves as the order parameter to quantify the first transition from the low-Wi symmetric state to either the D or C states (Fig. 2b,e). $I^{\prime}=0$ when the average flow through the upper and lower gaps equals the flow through the center. Transition to the D state results in $I^{\prime}>0$, since $\bar{u}_{0}$ decreases. Transition to the C state results in $I^{\prime}<0$, since $\bar{u}_{0}$ increases. $I^{\prime\prime}$ serves as the order parameter to quantify the second transition between the D and the AD states. $I^{\prime\prime}=0$ in the D state, since $\bar{u}_{+}=\bar{u}_{-}$ (Fig. 2b). In the AD state, fluid flows preferentially through either the upper or lower gap (Fig. 2c₁,c₂), resulting in $I^{\prime\prime}>0$ or $I^{\prime\prime}<0$, respectively.

The asymmetry parameters $I^{\prime}$, $I^{\prime\prime}$ are shown vs Wi in Fig. 3 for various values of $G$ and are fitted with a quartic (double-well) Landau-type potential minimized as:

where $\textnormal{Wi}_{c}=\textnormal{Wi}_{1}$ or $\textnormal{Wi}_{2}$ is the critical Weissenberg number for the bifurcation, and the order parameter $\epsilon$ can be either $I^{\prime}$ or $I^{\prime\prime}$. In all the fits, the growth rate coefficient $g$ is order unity, and the asymmetric term in $h$ quantifies system imperfections that bias a transition to a favored branch. The phenomenological Landau model for equilibrium phase transitions has long been found to provide a good description of bifurcation phenomena in driven nonequilibrium systems including Newtonian and viscoelastic flows [37, 38, 4, 6]. Eq. 2 describes forward (supercritical) pitchfork bifurcations without hysteresis.

For small $G=0.500$, the first transition in $I^{\prime}$ (Fig. 3a) occurs at $\textnormal{Wi}_{1}\approx 13$, and is a slightly imperfect ($h\approx-0.016$) supercritical pitchfork bifurcation where the favored branch ($I^{\prime}>0$) gives diverging (D) flow. The unfavored ($I^{\prime}<0$) branch was never observed, but its hypothetical existence is indicated in Fig. 3a by the dotted line. With increasing Wi, $I^{\prime}\rightarrow 1$, implying that almost no fluid passes between the two cylinders (as qualitatively evident from Fig. 2b). The second transition in $I^{\prime\prime}$ (Fig. 3b) from the D state to the asymmetric diverging (AD) state, occurs at $\textnormal{Wi}_{2}\approx$ 44. The imperfection in this second bifurcation is very small ($h\approx-0.0013$). The favored branch gives $I^{\prime\prime}>0$, but the unfavored $I^{\prime\prime}<0$ branch can also be reached and followed by initiating the flow at a high Wi and subsequent quasistatic reduction. The complete bifurcation diagram showing the total asymmetry $I=I^{\prime}+I^{\prime\prime}$ vs Wi for $G=0.500$ is shown in Fig. 3c. The first bifurcation results in $I\rightarrow 1$. The second bifurcation splits $I$ into two branches, $I\rightarrow 1.5$ or $I\rightarrow 0.5$.

The behavior for $G=0.603$ is shown in Fig. 3d-f. In this case, the first bifurcation occurs at $\textnormal{Wi}_{1}\approx 20$ (Fig. 3d). The asymmetric term in Eq. 2 is positive ($h\approx 0.033$), resulting in a preferred transition from symmetric to converging (C) flow. With increasing Wi, $I^{\prime}\rightarrow-1$, indicating that nearly all of the fluid passes between the cylinders (see Fig. 2e,f). Since $h$ is relatively large, the negative $I^{\prime}$ branch is strongly preferred. The positive $I^{\prime}$ branch (dotted line in Fig. 3d) is never observed experimentally. When the system selects the C state at the first transition, a second bifurcation is not observed, and $I^{\prime\prime}\approx 0$ for all Wi (Fig. 3e). The complete bifurcation diagram for $G=0.603$ is shown in Fig. 3f: since $I^{\prime\prime}\approx 0$, $I\approx I^{\prime}$.

The data shown in Figs. 2 and 3a-f demonstrate two disparate flow behaviors that are sensitive to the value of $G$. The first bifurcation to either the D or C states is well described as a supercritical pitchfork bifurcation quantified by $I^{\prime}$. The two states are different branches of the same bifurcation and the value of $G$ determines which branch is selected by changing the sign of the asymmetric term ($h$) in Eq. 2. This implies the existence of a specific intermediate value of $G$ at which the bifurcation of $I^{\prime}$ should be perfect ($h=0$) and the D or C states are equally likely.

By examining a range of intermediate values $0.56<G<0.60$, we confirmed this assumption, as exemplified by Fig. 3g-i for $G=0.588$. Here, increasing Wi quasistatically from 0, the favored positive $I^{\prime}$ branch (D state) is observed (Fig. 3g, triangles) on exceeding $\textnormal{Wi}_{1}\approx 20$. However, the imperfection is sufficiently small ($h\approx-0.007$), that by quasistatic reduction of Wi from a high value, the unfavored C state branch (squares) can also be followed. From the D state, on exceeding $\textnormal{Wi}_{2}\approx 66$, the second bifurcation to the AD state occurs to either positive or negative $I^{\prime\prime}$ with almost equal likelihood in a given experiment (Fig. 3h triangles) since $h\approx 0$. The complete bifurcation diagram for $G=0.588$ in Fig. 3i shows the bistable coexistence of the C and D states for $\textnormal{Wi}_{1}<\textnormal{Wi}<\textnormal{Wi}_{2}$. For $\textnormal{Wi}>\textnormal{Wi}_{2}$, the system is tristable, where the two AD branches and the C branch coexist.

The behavior observed in all eleven microchannels is summarized in a flow stability diagram in Wi–$G$ state space (Fig. 4). For channels with ‘small’ $G\lesssim 0.560$ the systems behave as exemplified by Figs. 2a-c and 3a-c, with a first bifurcation at $\textnormal{Wi}_{1}$ that is biased to the D state and a second bifurcation at $\textnormal{Wi}_{2}$ to the bistable AD state with either positive or negative $I^{\prime\prime}$. For ‘large’ $G\gtrsim 0.595$ the systems behave as shown by Figs. 2d-f and 3d-f; the bifurcation at $\textnormal{Wi}_{1}$ is biased to the C state (which is maintained until the onset of time-dependence at $\textnormal{Wi}>>\textnormal{Wi}_{1}$). For a narrow range of $0.560\lesssim G\lesssim 0.595$, beyond $\textnormal{Wi}_{1}$ it is possible to pass through a region of bistability between the D and C states before entering a tristable region beyond $\textnormal{Wi}_{2}$ comprising the bistable AD state and the C state. It is assumed that the D/C bistable region meets the low-Wi symmetric region at a hypothetical single point ($G_{c},\textnormal{Wi}_{1}$) in the state space (red crossed circle, Fig. 4) where the first bifurcation would be perfect ($h=0$). To the left of this point, the diagonal boundary between the D/C and D regions reflects the increasing imperfection of the first bifurcation with decreasing $G$. However, based on our experiments, to the right of $G_{c}$ the boundary between the D/C and the C states appears to be extremely abrupt. A small increase in $G>G_{c}$ causes a significant bias to the C state in the first bifurcation. A special case arises for $G=0.565$, where the tristable AD/C region is reached by passing through the bistable AD region, avoiding the bistable D/C region. Apparently, at this value of $G$ and Wi, the selected flow path can spontaneously switch from an edge to the center gap, and vice-versa for decreasing Wi.

Our extensive microfluidic experiments reveal the complex dynamical behavior of viscoelastic flow past side-by-side microcylinders with variable spacing. We observe flow transitions that are well-described by the Landau model as supercritical pitchfork bifurcations. If the intercylinder spacing is small, a first bifurcation is biased to favour a diverging flow state where the fluid avoids the intercylinder gap. From here, a second bifurcation leads to an asymmetric-divergent state where flow through just one of the two side gaps is preferred. The overall behavior is reminiscent of shear-thinning viscoelastic flow around a single obstacle [6, 7]. For large intercylinder spacing, the system undergoes a single bifurcation that is biased to a converging flow state where all of the fluid passes between the cylinders. For a small range of intermediate intercylinder gaps, it is possible to identify a region of bistability where the converging and diverging states coexist, which becomes a tristable region when the diverging state undergoes the second bifurcation.

This work is the first study of the creeping viscoelastic flow past side-by-side cylinders. Despite the modest increase in geometrical complexity from flow around a single cylinder, the dynamical behavior is significantly richer and provides an important and neglected stepping stone towards understanding viscoelastic flows through more complex arrays of objects. Our results suggest a new interpretation for how viscoelastic fluids select preferred flow paths through ordered and disordered porous arrays (e.g., Ref. 15), indicating that each individual obstacle should be considered as a bifurcation point, which can be perfect (ordered) or imperfect (disordered), depending on the spacing between array elements.

We gratefully acknowledge the support of the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan, and also funding from the Japan Society for the Promotion of Science (JSPS, Grant Nos. 18K03958, 18H01135 and 20K14656) and the Joint Research Projects (JRPs) supported by the JSPS and the Swiss National Science Foundation (SNSF).