Elastic instability of wormlike micelle solution flow in serpentine channels

Abstract

Wormlike micelle (WLM) solutions are abundant in energy, environmental, and industrial applications, which often rely on their flow through tortuous channels. How does the interplay between fluid rheology and channel geometry influence the flow behavior? Here, we address this question by experimentally visualizing and quantifying the flow of a semi-dilute WLM solution in millifluidic serpentine channels. At low flow rates, the base flow is steady and laminar, with strong asymmetry and wall slip. When the flow rate exceeds a critical threshold, the flow exhibits an elastic instability, producing spatially-heterogeneous, unsteady three-dimensional (3D) flow characterized by two notable features: (i) the formation and persistence of stagnant but strongly-fluctuating and multistable “dead zones” in channel bends, and (ii) intermittent 3D “twists” throughout the bulk flow. The geometry of these dead zones and twisting events can be rationalized by considering the minimization of local streamline curvature to reduce flow-generated elastic stresses. Altogether, our results shed new light into how the interplay between solution rheology and tortuous boundary geometry influences WLM flow behavior, with implications for predicting and controlling WLM flows in a broad range of complex environments.

---

1. Introduction

Above a threshold concentration in solution, many surfactants self-assemble into long, flexible, polymer-like chains reaching microns in length. However, unlike polymers, these wormlike micelles (WLMs) constantly break and reform at thermal equilibrium and under flow. Hence, they are often termed “living polymers”.^1,2 When concentrated enough to be semi-dilute, WLM solutions are both strongly shear thinning and highly elastic, and therefore, are used in a wide range of energy and industrial applications, such as subsurface proppant transport,^3–6 groundwater remediation,^7,8 heating and cooling,^9–11 drag reduction,^12–16 lubrication,¹⁷ drug manufacturing and delivery,^18–20 and as rheological modifiers in consumer products.^4,14,21,22

All of these applications involve the flow of WLM solutions through complex, confined spaces, such as winding pipes and porous rocks, where the boundary geometry introduces tortuosity to the flow. This tortuous flow in turn couples to the complex WLM solution rheology to produce unusual flow behaviors. For example, in Couette-type geometries, the homogeneous base flow can become unstable, separating into low and high shear bands that together maintain a constant shear stress with increasing shear rate^23–28—a form of shear localization. The high shear rate band can subsequently exhibit another instability that produces unsteady flow,^26,29–33 with pronounced spatiotemporal fluctuations similar to those exhibited by elastic polymer solutions in tortuous spaces.^34–40 This instability arises at low Reynolds number (Re ≪ 1)—so the effects of inertia, which often causes flow instabilities, are negligible. Instead, it is generated by fluid elasticity, and is therefore often parameterized by the Weissenberg number, Wi, which compares the relative strength of elastic and viscous stresses. This purely-elastic instability has also been documented for WLM solutions flowing in a range of microfluidic geometries, such as through rectangular channels,^41,42 around sharp bends,^43,44 through contractions and expansions,^45,46 in cavities,⁴⁷ and past a single obstacle^48–55 or two-dimensional (2D) array of obstacles.^56–58 However, such geometries have multiple confounding factors, such as mixed shear and extensional flow topology and the presence of stagnation points at solid surfaces, that can simultaneously influence the flow behavior. Thus, a clear understanding of how the coupling between tortuous boundary geometry and fluid rheology influences the complex flow behavior of WLM solutions is still lacking.

Here, we simplify this problem by studying WLM solution flow at Re ≪ 1 and Wi ∼ 1 in serpentine channels comprising successive semi-circular half-loops. This microfluidic geometry enables us to isolate the role of fluid streamline curvature—which is known to be a key factor controlling purely-elastic flow instabilities since it generates elastic hoop stresses^36,37—in shaping the unstable flow. Indeed, our work is inspired by prior studies of elastic instabilities in polymer solutions;^59–64 these provide useful intuition, but are not directly reflective of the flow behavior of WLM solutions, which have fundamentally different microstructures and rheological properties. We find that the base steady, laminar flow exhibits strong asymmetry and wall slip at the lowest flow rates tested due to the unique rheology of the WLM solution, which enables it to support shear localization. Above a critical Wi = Wi_c, the flow becomes unstable, producing spatially-heterogeneous unsteady three-dimensional (3D) flow. This unsteady flow has two notable features: (i) persistent, fluctuating dead zones, which form due to the ability of the WLM solution to support shear localization, and (ii) slow “twisting” events in 3D. Our characterization of these unusual flow behaviors suggests that when WLM solutions are forced to flow in confined, tortuous spaces, they develop pathways of least resistance that minimize local streamline curvature. Altogether, our work provides a key first step toward making sense of complex fluid flows in complex spaces, potentially enabling their prediction and control in a broad range of applications.

2. Methods

2.1. Test fluid and preparation

We use a well-characterized semi-dilute, entangled WLM solution composed of 100 mM cetylpyridinium chloride (CPyCl; Sigma-Aldrich) and 60 mM sodium salicylate (NaSal; Sigma-Aldrich) in ultrapure Millipore water.^{45–47,51,52,58,65–70} We prepare the solution by first dissolving the CPyCl, which is a cationic surfactant, and then adding the NaSal, an aromatic organic salt, and stirring the solution for 24 h. The addition of NaSal induces the elongation of micelles due to electrostatic screening of the charged headgroups and steric interaction with the aromatic ring of the organic salt.⁷⁰ We let the resulting WLM solution rest for >48 h to allow air bubbles to disappear before use. We also add 1 μm diameter fluorescent, amine-functionalized polystyrene tracer particles (Invitrogen; yellow-green) to the solution at 20 ppm for flow visualization and particle image velocimetry (PIV).

2.2. Shear rheology

To characterize the shear rheology of the WLM solution, we use a stress-controlled Anton Paar MCR501 rheometer fitted with a truncated cone-plate geometry (CP50-2: 50 mm diameter, 2°, 53 μm) at a controlled temperature of 25 °C. We measure steady-state flow curves at fixed values of the shear rate [small gamma, Greek, dot above], ramping up and down between 0.01 and 10 s⁻¹; shear rates exceeding 10 s⁻¹ result in fluid ejection due to an elastic instability. The results, averaged over three replicate measurements, are shown in Fig. 1(a) and (b). As shown in Fig. 1(a), the shear stress σ increases approximately linearly with shear rate up to [small gamma, Greek, dot above]₀ ≈ 0.3 s⁻¹, and then plateaus at a value σ* = 13 Pa, a signature that may reflect the presence of shear banding.^24–27 As we show later on, this strong shear thinning of the WLM solution has an important consequence for its flow in serpentine channels: it enables them to support shear localization, i.e., separate into low and high shear regions under a constant imposed bulk flow rate. As shown in Fig. 1(b), the corresponding dynamic viscosity η([small gamma, Greek, dot above]) = σ([small gamma, Greek, dot above])/[small gamma, Greek, dot above] is nearly constant at η₀ = 49 Pa s at low shear rates, and then continually decreases (shear thinning) with increasing [small gamma, Greek, dot above] > [small gamma, Greek, dot above]₀. This shear thinning behavior is described well by the Carreau–Yasuda model,⁷¹ as demonstrated in previous studies of similar or identical WLM solutions;^{49,51,54,55,58,69,72} the model is shown by the solid lines: η([small gamma, Greek, dot above]) = η_∞ + (η₀ − η_∞)(1 + ([small gamma, Greek, dot above]/[small gamma, Greek, dot above]₀)^a)^(n−1)/a, where a = 3.8 describes the transition to the shear thinning regime, the power law index n ≈ 0, and we set η_∞ equal to the solvent viscosity of 1 mPa s given the absence of a measurable high shear rate Newtonian plateau. Importantly, the Carreau–Yasuda model does not capture any variations in local fluid viscosity due to any shear banding behavior; we simply use it to identify the onset and strength of shear thinning in our flow experiments.

Fig. 1 Shear rheology of 100/60 mM CPyCl/NaSal WLM solution. (a) Shear stress σ and (b) shear viscosity η as a function of imposed shear rate [small gamma, Greek, dot above], exhibiting a shear stress plateau. Solid lines correspond to a Carreau–Yasuda fit. (c) First normal stress difference N₁ as a function of imposed shear rate [small gamma, Greek, dot above], with power law fit (solid line). (d) Storage G′ and loss G′′ moduli as a function of angular frequency ω under small amplitude oscillatory shear in the linear viscoelastic regime at 1% strain amplitude. Solid curves are fits to the single-mode Maxwell model with a characteristic relaxation time λ_M = 1.4 s. Error bars correspond to the standard deviation over 3 replicate measurements.

To characterize the solution elasticity, we also measure the first normal stress difference N₁ and linear elastic storage and loss moduli G′ and G′′, respectively. As shown in Fig. 1(c), N₁ increases with shear rate as a power law (solid curve), , with and power law index n_N1 = 0.85. The instrument resolution due to normal force sensitivity is roughly N₁ ≈ 10 Pa. Thus, estimates of N₁ using the power law model at shear rates below [small gamma, Greek, dot above] ≈ 1 s⁻¹ are extrapolations of this power law and assume that the power law dependence holds in the non-shear thinning regime at low shear rates. Given that the measured power law dependence of N₁ occurs over a range of shear rates corresponding to the shear stress plateau in Fig. 1(a), an interesting direction for future work could be to explore how this power law dependence relates to, and potentially can help to estimate the amount of, shear banding. As shown by the curves in Fig. 1(d), G′ and G′′ determined from small-amplitude oscillation over a range of angular frequencies ω are well described by a single-mode Maxwell model: and , where G₀ = 29 Pa is the plateau modulus and λ_M = 1.4 s is the Maxwell relaxation time. We thereby estimate the mesh size of the entangled micellar network as .⁷³ Based on prior estimates of the persistence length for the same WLM solution (l_p ≈ 28 nm^67,70), we estimate the average micelle contour length to be ∼830 nm.^70,74 The rheology measurements also enable estimates of the micelle breakage and reptation time scales τ_break = 0.49 s and τ_rept = 3.5 s, respectively, using the framework of Turner and Cates.⁷⁵ Their ratio , confirming that the micelles are in the fast-breaking limit where the fluid is well-described by the idealized Maxwell model.⁷⁵

2.3. Device architecture and fabrication

The serpentine channels used for the flow experiments comprise 19 semi-circular half-loops with inner radius R_i, fixed channel width W = 1 mm, outer radius R_o = R_i + W, fixed channel height H = 2 mm, and fixed aspect ratio H/W = 2; an optical top-down image of a representative example is shown in Fig. 2(a). We choose a large channel aspect ratio to ensure that the dominant velocity gradients arise in the flow imaging plane. All the quantitative results presented in this paper are for channels with R_i = 300 μm; however, as shown in ESI† Movies S3 and S4, we observe similar behavior in channels with R_i = 500 and 1000 μm.

Fig. 2 Serpentine channel geometry and experimental setup. (a) Brightfield optical image of the serpentine channel device comprised of 19 semi-circular half-loops with inner radius R_i and channel width W. (b) Schematic of experimental setup for flow visualization.

We design the body of each channel [dark blue in Fig. 2(b)] using CAD software (Onshape) and 3D-print it with a clear methacrylate-based polymeric resin (FLGPCL04) using a FormLabs Form 3 stereolithography printer. This resin is inert and likely does not interact with or alter the WLM solution through specific chemical interactions.^2,76,77 Moreover, while slight surfactant adsorption to the channel walls may occur, given the device geometry, we estimate that the resulting change in the solution concentration is at most 0.05%—negligible compared to the concentration used in our experiments. We assemble each millifluidic device by screwing together this 3D-printed channel body and a laser-cut transparent acrylic top sheet (McMaster-Carr) [gray in Fig. 2(b)], with a ∼1 mm-thick layer of polydimethylsiloxane (PDMS; Dow SYLGARD 184, 8.5 : 1.5 base : curing agent by weight) sandwiched in between to act as a gasket. We then glue flexible Tygon tubing (McMaster-Carr) in the inlet and outlet using a 2-part watertight epoxy (JB MarineWeld).

2.4. Flow imaging

Before each flow experiment, we flush the device to be used with ultrapure Millipore water to remove any residual debris or air bubbles. We then fully saturate the channel with the WLM solution to be tested at a flow rate of 1 mL h⁻¹ for at least 3 h using a syringe pump (Harvard Apparatus PHD 2000). The device is affixed to the stage of an inverted Nikon A1R+ laser scanning confocal fluorescence microscope, ensuring that the syringe pump, device, and outlet waste jar are at the same height to avoid hydrostatic pressure differences, as shown in Fig. 2(b).

In each experiment, we impose a stepwise ramp of increasing flow rate Q from 0.5 to 6 mL h⁻¹, allowing the flow to equilibrate at each flow rate for 30 min before commencing imaging. For each flow rate tested, we use a 4× objective lens to acquire 2 min long (∼86 λ_M) movies at 30 frames per second focused on a horizontal optical plane, 37 μm in thickness, centered at the mid-channel height to avoid boundary effects induced by the top and bottom channel walls. Each movie is taken at each successive half-loop along the length of the serpentine array, with a 3167 × 3167 μm² field of view at a resolution of 6 μm per pixel. The fluorescent tracer particles seeded in the fluid are excited by a 488 nm laser and their emission is detected using a 500–550 nm sensor.

To characterize the flow behavior, we define a nominal shear rate using the channel half-width W/2 as the characteristic length scale: . Our experiments probe the range 0.14 ≤ [small gamma, Greek, dot above] ≤ 1.7 s⁻¹; as shown by the gray region in Fig. 1(a), this range corresponds to strongly shear thinning conditions for the WLM solution. We thereby define the Weissenberg number Wi, which compares the relative strength of elastic and viscous stresses, as , where N₁([small gamma, Greek, dot above]) and η([small gamma, Greek, dot above]) are given by the power law and Carreau–Yasuda fits to the bulk shear rheology, respectively. Our experiments probe the range 0.17 < Wi < 0.78, and the onset of shear thinning arises at Wi₀ = Wi([small gamma, Greek, dot above]₀) = 0.19. For flow conditions exceeding Wi = 0.78, we are unable to obtain PIV measurements for quantitative analysis.

To visualize the WLM solution flow, we generate fluid pathline images and movies using grouped projections of mean pixel intensity over successive frames using the open-source software ImageJ.^78,79 Pathlines correspond to trajectories of tracer particles as they are advected by the fluid flow. Example movies are shown in ESI† Movies S1–S4. To quantitatively analyze the flow field, we directly measure each time-resolved 2D velocity field u using PIVlab,⁸⁰ manually isolating regions in the imaged field of view that only contain fluid flow. We then use the FFT window deformation algorithm with 4 passes of decreasing interrogation area from 64–32–16–16 pixel box sizes and 50% overlap, producing a spatial grid resolution of 50 μm for which particles in each PIV voxel travel approximately 4–8 pixels between successive frames.^81,82 Given that the channel width of 1 mm greatly exceeds the spatial PIV resolution, this approach resolves flow gradients in the selected 2D imaging plane. Similarly, the imaging frame rate (30 frames per second) sets the temporal resolution of our velocity fields (33 ms). Given the fluid relaxation time of 1.4 s and total imaging time at each location of 120 s, our approach enables us to resolve temporal variations in the flow field for variations occurring over a broad range of timescales from 33 ms to 120 s. Finally, we use in-house MATLAB codes to process the velocity fields for all analyses. We characterize the flow in each channel using three metrics: the normalized velocity magnitude |u|/max(|u|), root mean square velocity fluctuations and shear rate , where x and t represent the 2D position in the horizontal optical plane and time, respectively, and is the rate-of-strain tensor computed from the measured velocity field.

3. Results and discussion

3.1. Base steady laminar flow

We first characterize the steady base flow of the WLM solution in the serpentine channels at low Wi. A representative pathline image of the laminar steady flow at Wi = 0.18 is shown in Fig. 3(a), and the velocity profiles at the apex and exit of a half-loop are shown for the regions outlined in yellow and red in Fig. 3(b), respectively; the corresponding ESI† Movie S1 shows that the pathlines remain unchanging over time. Notably, at the apex of the bend, the velocity profile is highly asymmetric and skewed toward the outer bend. While the curvature of the channel walls is known to induce slight asymmetry in Newtonian flows,^61,63,64,83 in this case the velocity profiles are starkly different. In particular, as shown in Fig. 3(b), the flow is plug-like, unlike the typical parabolic profile associated with Poiseuille flow, because of the shear thinning of the WLM solution. Additionally, near the channel walls, the velocity does not decrease to zero—indicating strong slip caused by shear localization at the fluid–solid boundary, as has been reported for WLM flows in other geometries.^17,26,66,69 We estimate characteristic slip velocities by linearly extrapolating the velocity profiles near the solid surfaces, finding average slip velocities of 0.07, 0.22, and 0.48 mm s⁻¹ at flow rates Q = 0.5, 1.5, and 3 mL h⁻¹, respectively, indicating that the extent of slip increases monotonically with increasing flow rate in the steady flow regime. The corresponding slip lengths (>100 μm in all cases) greatly exceed the spatial resolution of PIV (50 μm), which additionally confirms the presence of strong slip. The estimated slip velocities are persistent in time in the steady flow regime, indicating sustained shear localization occurring at the channel walls.

Fig. 3 Base steady laminar flow exhibits strong shear thinning and slip at surfaces. (a) Pathline image of flow in the R_i = 300 μm device at Wi = 0.18 < Wi_c, generated using a moving average of fluorescence intensity over 60 frames. Scale bar is 500 μm. (b) Time-averaged (over 86λ) velocity profiles, normalized by the maximum velocity, at the apex (yellow rectangle) and end (red rectangle) of a half-loop. All velocity profiles exhibit highly asymmetric plug flow and indicate the presence of a substantial slip velocity at surfaces. The PIV grid resolution used is 0.05 W. (c) Streamline visualization and color maps of the normalized velocity magnitudes at selected half-loop locations further reflect strong plug flow and asymmetry in the base laminar flow state.

Interestingly, we observe hints of shear inhomogeneity reminiscent of shear banding behavior—referred to here as “shear localization”—when Wi is increased slightly to 0.25, as shown by the inflection in the blue profile shown in the right-hand plot of Fig. 3(b). This effect diminishes upon further increasing Wi, as shown by the green profile in Fig. 3(b), suggesting that the flow re-homogenizes across the channel width. A representative image showing the velocity magnitudes and flow streamlines across multiple half-loops is shown in Fig. 3(c) for this case of Wi = 0.43, at which the flow remains laminar and steady. The analyses in Fig. 3 present 2D velocity profiles measured at the channel mid-plane; however, due to the 3D nature of the channel, steady secondary flows not captured by our 2D measurement likely arise due to the inversion of curvature in the serpentine geometry, as have been reported previously.^83–86

3.2. Onset of 3D unsteady flow

Further increasing Wi gives rise to an elastic flow instability where this base steady flow becomes unstable to 3D unsteady flow. An example is shown in Fig. 4(a); for this channel geometry, the onset of the instability—determined as when the rms velocity fluctuations exceed a threshold value of 0.2 set by the noise threshold of the PIV measurements—occurs at a critical value Wi = Wi_c ≈ 0.7. Such purely-elastic instabilities are typically generated through the coupling of elastic normal stresses and streamline curvature, which gives rise to destabilizing hoop stresses.^34,36–38 This coupling can also be described using a dimensionless parameter developed by Pakdel and McKinley:^36,37, where U = Q/(HW) is the average flow velocity. In our experiments, the critical value of M = M_c ≈ 2.5, in good agreement with prior measurements across a broad range of fluids and geometries.³⁸

Fig. 4 Unsteady 3D flow arising from an elastic instability. (a) Color maps of the flow field measured at the channel mid-height in the unsteady flow state showing the (top) normalized velocity magnitude |u|/max(|u|), (middle) root-mean-square (rms) velocity magnitude fluctuations and (bottom) shear rate [small gamma, Greek, dot above](x) in the R_i = 300 μm channel at Wi = 0.78 > Wi_c. The top and bottom rows show near-instantaneous flow fields that are time-averaged over a duration of 1 λ_M to reduce noise. The middle row shows rms velocity fluctuations over a duration of 86 λ_M, highlighting the spatial heterogeneity and temporal persistence of flow fluctuations. All color maps comprise stitched fields-of-view obtained sequentially during flow imaging, giving rise to the visible discontinuities between adjacent half-loops. (b) Power spectral densities of the magnitude of velocity fluctuations indicate that above the onset of the elastic instability, more energy is dissipated in fluctuations across a broad range of frequencies, with no characteristic time scale.

As shown by the maps in Fig. 4(a), all three flow metrics exhibit pronounced heterogeneity throughout the channel. To further characterize the flow dynamics, we examine the temporal power spectral density (PSD) of the magnitude of velocity fluctuations across different frequencies f: , where L is the duration of the signal, F_s is the sampling rate, and FFT denotes a fast Fourier transform. Our results, spatially averaged over all PIV voxels, are shown in Fig. 4(b). Above the onset of the elastic instability, the PSD exhibits a broad power law decay E(f) ∼ f^−α, with α ≈ 2.6 at Wi = 0.78, reflecting energy dissipation by flow fluctuations across a broad range of time scales—consistent with previous characterizations of polymeric and WLM elastic instabilities in other confined geometries.^{26,30,39,45,46,63,64,87–89} The absence of distinct peaks in the PSD indicates that the unsteady flow is aperiodic with no characteristic time scales. Analysis of the instantaneous velocity profiles shows the presence of pronounced slip behavior consistent with the strong slip present in the steady flow regime (Section 3.1); however, the estimated slip velocities now vary greatly over time (e.g. ranging from 0 to 2.54 mm s⁻¹ at Q = 6 mL h⁻¹) due to the strongly fluctuating nature of the flow.

3.2.1. Dead zone formation, dynamics, and size are shaped by fluid rheology. Close inspection of the unstable flow shown in Fig. 4(a) reveals a fascinating consequence of the WLM solution's ability to support shear localization: As shown by the dark blue regions in the top row of Fig. 4(a), non-flowing “dead zones” (identified by |u|/max(|u|) < 0.05) form and persist in the downstream portion of some of the half-loops. The corresponding ESI† Movie S2 shows these dynamics using pathline imaging. While they are reminiscent of corner eddies reported for viscoelastic WLM solution flows around sharp bends,^43,44 the dead zones formed in our serpentine channels are truly near-stagnant regions with minimal flow, whereas eddies are characterized by non-zero circulating flow. This feature is illustrated by the shear rate profiles across the channel width (inset) in Fig. 5(a), and an example instantaneous pathline image with two dead zones demarcated by the yellow dashed outlines is shown in Fig. 5(b). As shown in the dark-colored profiles in Fig. 5(a), a clear peak of high shear rate separates the dead zone ([small gamma, Greek, dot above] ≈ 0 s⁻¹) and the flowing regions, generating a strong shear rate gradient at the dead zone boundary; by contrast, if the fluid could not support shear localization, dead zones would not form, as previously documented for non-shear banding elastic polymer solutions.^61–64

Fig. 5 Characterization of dead zone formation, dynamics, and size. (a) Temporal evolution of shear rate profiles along the channel width (location shown in inset). Successive curves are plotted at time intervals of 6.7 s. At early times (0 < t < 90 s), a large peak in shear rate separates the dead zone ([small gamma, Greek, dot above] ≈ 0 s⁻¹) from the bulk flow, while at later times the dead zone is washed away and the shear rate profile re-homogenizes. (b) Example instantaneous pathline image of the unsteady flow state in a selected half-loop. Yellow dashed lines demarcate example dead zone regions. Scale bar shows 500 μm. (c) Temporal variation of dead zone size A_DZ in selected half-loops shows considerable fluctuations. The shaded gray region denotes the range of theoretically-predicted dead zone areas from considering minimization of local streamline curvature. (d) Power spectral densities of dead zone size fluctuations indicate no characteristic time scale. (e) Probability density distribution of normalized dead zone areas A_DZ/A_half loop with increasing advective Deborah number De_adv shows a broad increase in dead zone sizes for De_adv > 1. (f) Degree of multistability ξ describing the range of observed dead zone sizes at a given flow condition shows onset of multistable dead zone behavior above De_adv = 1. All data in (a)–(d) are for Wi = 0.78.

To characterize their dynamics, we track variations in the dead zone areas (A_DZ) over time. Examples for selected half loops along the length of the channel are shown by the traces in Fig. 5(c). Notably, the dead zone sizes fluctuate dramatically over time, often persisting over durations much longer than the characteristic relaxation time λ_M = 1.4 s. These fluctuations do not have any clear periodicity; instead, they mirror the spectrum of bulk velocity fluctuations shown in Fig. 4(b), as reflected by the broad decay in the power spectral density E_DZ(f) in Fig. 5(d)—indicating that dead zone fluctuations are established by the unstable fluctuations in the freely-flowing fluid outside the dead zone. Despite these fluctuations, however, all dead zones are bounded by a maximal size A^max_DZ, indicated by the gray region in Fig. 5(c).

Intriguingly, for a given unstable flow rate and at a given time, not all half-loops contain dead zones—as shown in Fig. 4(a) and 5(c). Instead, they exhibit multistable behavior: some half-loops do not have dead zones (the dead zone-free state), while other half-loops have persistent dead zones of varying sizes (the dead zone-containing state). As time progresses, a given half-loop can randomly switch from the dead zone-free state to the dead zone-containing state and vice versa; an example is shown by the different traces in Fig. 5(a), in which a dead zone eventually washes away, as well as those in Fig. 5(c). This behavior is reminiscent of the multistability exhibited by eddies during the unstable flow of elastic polymer solutions in pore constriction arrays,^90–92 in which multistability arises when flow-induced elastic stresses do not have sufficient time to relax as they are advected between constrictions. Hence, following this prior work, we parameterize the onset of the multistability observed for WLM solutions in serpentine channels using a streamwise Deborah number comparing the stress relaxation time of the WLM solution to the time required for the fluid to be advected between adjacent half-loops: where t_adv = HA_half loop/Q and A_half loop is the half-loop area. Consistent with this prior work, we also find that multistability arises when De_adv ≳ 1. For example, Fig. 5(e) shows a color map of the probability density function (p.d.f.) of the dead zone areas aggregated over the different half-loops, normalized by A_half loop, for the different De_adv tested—showing that this p.d.f. abruptly begins to broaden, a hallmark of multistability, when De_adv ≳ 1. Furthermore, also following prior work, we define the degree of multistability , where ξ = 0 indicates no multistable behavior and ξ = 1 describes the maximum possible extent of multistable behavior. As shown in Fig. 5(f), we again see an abrupt increase in the degree of multistability ξ when De_adv ≳ 1. Taken together, these results demonstrate that fluid elasticity plays a key role in determining dead zone formation.

Motivated by these findings, as well as previous work describing the shapes of eddies formed in elastic polymer solutions entering contractions,^93–97 we conjecture that considering the elastic stresses generated during WLM solution flow in a serpentine channel can help predict the maximal dead zone size A^max_DZ. In particular, given that these elastic hoop stresses are generated by streamline curvature,^36,37 we expect that the interface between a dead zone and the surrounding freely-flowing region forms to minimize local streamline curvature. This idea can be formulated using a simple geometric description (detailed in the Appendix): for a given half-loop, we expect that the largest dead zone that can form in a half-loop is bounded by a straight line tangent to the inner bend of the loop, as shown in Fig. 7. This prediction, shown by the gray shaded region in Fig. 5(c), agrees reasonably well with the experimental measurements as an upper bound A^max_DZ to the individual traces of A_DZ.

Fig. 6 3D twisting flows intermittently arise in the unsteady flow. (a) Snapshots of the velocity field (arrows) superimposed on a color map of normalized velocity magnitude at different times, showing a representative twisting event. (b) Probability density function of the divergence of the measured 2D flow field suggests an increase in 3D flow above the onset of instability (Wi > Wi_c ≈ 0.7). (c) Twisting reduces the hydrodynamic tortuosity of fluid streamlines: the blue dashed curve is the center semi-circular streamline in the base flow, and the green dashed curve is the projected 2D sinusoidal path enabled by 3D twisting. Pathline images are generated using a moving average of fluorescence intensity over 10 frames. Scale bar is 500 μm. Data in (a) and (c) are for Wi = 0.78.

Fig. 7 Schematic describing the geometric prediction of dead zone size and shape in the serpentine channels.

3.2.2. 3D twisting flow events reduce local hydrodynamic tortuosity. ESI† Movie S2 also reveals another unusual feature of the unstable WLM solution flow: 3D flow inversion events characterized by lower velocity “twists” in the imaging plane. Additional examples are shown in ESI† Movie S4. Fig. 6(a) presents snapshots taken at different times of the fluid velocity field (arrows, with the normalized velocity magnitude shown by the color map) corresponding to one such twisting event. The initial faster-flowing, primarily 2D flow (bright green/yellow/red region in the first panel) is overtaken by a slower-moving flow (dark blue region) that twists in the z-direction, as reflected by a curvy “disclination” of the velocity vectors. Note that these vector fields correspond to a 2D representation of the flow in the selected imaging plane and do not show the out-of-plane velocity component. Fig. 6(b) confirms the presence of 3D flow in this unsteady flow state by plotting the probability density function of the divergence of the 2D velocity field, , using the in-plane velocity components. Under steady flow conditions (Wi = 0.18–0.43), the 2D divergence is close to zero, with slight deviations due to small weakly-3D secondary flows generated by the inversion of curvature at corners and bends;^83–86 however, in the unstable flow state (Wi = 0.78), an appreciable fraction of the 2D divergence becomes non-zero, reflecting flow in the third dimension.

Despite their complexity, we attempt to rationalize the basic geometric properties of these 3D inversion events by again considering minimization of the local streamline curvature—which here manifests as a reduction of the 2D hydrodynamic tortuosity, τ. In the base laminar flow, the streamlines can be approximated as a semi-circular path shown by the blue dashed line in Fig. 6(c); the associated tortuosity, given by the ratio of the total path length to the horizontal streamwise distance traveled, is then . Guided by the visualization in ESI† Movie S4 and Fig. 6(a), we approximate the “disclination” boundary of the twisting flow as a periodic sinusoidal function with amplitude R_i and period 2(R_i + R_o). As shown by the green dashed line in Fig. 6(c), this approximation is in good agreement with the experimental observations. The resulting 2D hydrodynamic tortuosity is then given by , where x and y represent the streamwise and transverse position coordinates. For this channel, this calculation then yields a tortuosity of 1.08, representing a 31% reduction from the laminar base case.

4. Conclusions

In this study, we investigated the flow behavior of a highly elastic, shear thinning, semi-dilute WLM solution in serpentine channels at low Reynolds number (Re ≪ 1) and moderate Weissenberg numbers (Wi ∼ 1). The serpentine microchannels are composed of successive semi-circular half-loops, thereby enabling us to experimentally isolate the influence of streamline curvature—a fundamental feature of many real-world flows—on the flow behavior without being confounded by other geometric complexities. Our flow visualization experiments revealed three key phenomena:

1. At low Wi, the base flow is steady and laminar but exhibits spatial asymmetry with wall slip, reflecting the shear thinning and shear localization tendencies of the WLM solution. Above a critical Wi = Wi_c ≈ 0.7, the flow undergoes an elastic instability and transitions to a 3D unsteady flow state characterized by pronounced spatiotemporal velocity fluctuations. This transition occurs at a critical Pakdel–McKinley parameter value of M_c ≈ 2.5, consistent with other measurements across different viscoelastic fluids and geometries.

2. Alongside this unstable bulk flow, dead zones of stagnant fluid form in the downstream portion of half-loops—reflecting the ability of the WLM solution to support shear localization, complementing reports of dead zone formation for other types of complex fluids.^98–102 Due to coupling to the velocity fluctuations in the bulk flow, these dead zones fluctuate in their size; however, they are bounded by a maximal size that minimizes the fluid streamline curvature, and therefore the generation of elastic stresses. Dead zones also exhibit multistable behavior—forming and persisting in some half-loops, not forming in other half-loops, and randomly switching between these two states. This multistability arises when the advective Deborah number De_adv ≳ 1, indicating that it occurs because elastic stresses generated by fluid flow do not have sufficient time to relax between consecutive half-loops. This finding expands on previous studies of multistability, which were restricted to elastic polymer solutions flowing through pore constriction arrays,^90–92 to a broader range of fluids and flow geometries.

3. The unstable flow state also features intermittent, 3D “twisting” velocity inversion events amid the spatiotemporally-fluctuating bulk flow. These twisting events reduce the hydrodynamic tortuosity compared to the base flow state, and their geometric structure can also be rationalized as minimizing the fluid streamline curvature, and therefore the generation of elastic stresses.

These findings highlight how the shear localization and elastic properties of WLM solutions give rise to unusual flow behaviors in tortuous channels, expanding current understanding of complex fluid flows in complex geometries. Further investigating the physics underlying the formation of the dead zones and 3D twisting events revealed by our experiments will be a useful direction for future work.

Several extensions of our study offer opportunities for future research. While our imaging focused on a single optical plane at the mid-channel height, full 3D velocimetry in channels of varying aspect ratio (H/W) would provide more comprehensive insights into the nature of the 3D twisting events we observed. Furthermore, direct measurements of the stress field during flow could provide deeper insights into the coupling between elastic stresses and flow structures (e.g., dead zones, 3D twisting events). Our experiments probed a wide range of shear rates, but they all fell on the stress plateau measured in the bulk shear rheology shown in Fig. 1(a); exploring a wider range of shear rates below and above this plateau will help further elucidate WLM flow dynamics across a broader range of conditions.

Finally, while here we focused on a single semi-dilute, linear, entangled WLM solution, future studies could explore other formulations with different rheological properties (e.g., varying the degrees of shear thinning and elasticity by tuning the molar ratio of surfactant to salt). Studies could also probe how variations in the micellar microstructure^{67,70,103,104} influence the flow behavior—which could potentially guide ways to use stimulus-responsive WLMs^10,18 to achieve targeted local flow behaviors. For example, controllably forming stagnant dead zones could be a mechanism to redirect flow on demand, while conversely, prevention of dead zone formation could help promote fluid homogenization.

Author contributions

E. Y. C.: conceptualization, methodology, investigation, formal analysis, writing; S. S. D.: conceptualization, resources, supervision, writing.

Data availability

The data that support the findings of this study are available within the article. Supplementary movies are available online at https://github.com/emchen6639/WLMserpentine.

Conflicts of interest

There are no conflicts to declare.

Appendix

Geometric prediction of dead zone shape and size

For one unit cell of the serpentine geometry, the upper boundary can be described as
where x and y represent the streamwise and transverse position coordinates. Similarly, the bottom boundary can be described as

Given a selected fixed point (x_b, y_b) on the outer boundary of a half loop, as depicted in Fig. 7, we seek to determine the location on the following inner bend (x*, y*) such that the line between the two points is tangent to the inner half-loop boundary, thereby minimizing local streamline curvature. The slope between the two points is , which is equal to the derivative of inner curve parameterization evaluated at (x*, y*): . Further, the point (x*, y*) must lie on both the inner curve equation and the tangent line equation. Accordingly, (x*, y*) is straightforwardly determined from the solution of:

Note that only one of two solutions is physical in the context of our experiments. By varying the selected point (x_b, y_b), we obtain predictions of the dead zone shape and size where the dead zone is bounded by the resulting tangent line and the outer channel boundary. Example solutions are shown by the red lines in Fig. 7. The maximum dead zone size is set by the physical constraints of the bounding geometry as conservation of mass requires that some finite region of the channel must contain non-zero flow velocity. The resulting range of theoretical dead zone size predictions 0 ≤ A^pred_DZ ≤ A^max_DZ is given by the shaded region in Fig. 5(c) and is in good agreement with the experimentally-measured values of dead zone sizes.

Acknowledgements

We acknowledge support from the Princeton Center for Complex Materials (PCCM), a National Science Foundation (NSF) Materials Research Science and Engineering Center (MRSEC; DMR-2011750), as well as the Camille Dreyfus Teacher-Scholar Program of the Camille and Henry Dreyfus Foundation. We also acknowledge the use of the Imaging and Analysis Center (IAC) operated by the Princeton Materials Institute at Princeton University. It is a pleasure to thank Howard Stone for useful feedback.

References

1. R. G. Larson, The Structure and Rheology of Complex Fluids, Oxford University Press, 1999.
2. C. A. Dreiss, Wormlike Micelles: Advances in Systems, Characterisation and Applications, The Royal Society of Chemistry, 2017, pp. 1–8.
3. B. Chase, W. Chmilowski, R. Marcinew, C. Mitchell and Y. Dang, Oilfield Rev., 1997, 9, 20–33.
4. R. Zana and E. W. Kaler, Giant Micelles Properties and Applications, CRC Press, Boca Raton, 1st edn, 2007.
5. E. Jafari Nodoushan, T. Yi, Y. J. Lee and N. Kim, Fluids, 2019, 4, 173.
6. A. Du, C. Ye, J. Chen, F. Sun and J. Mao, Colloids Surf., A, 2023, 672, 131707.
7. J.-H. Kwon, M.-K. Ji, R. Kumar, M. M. Islam, M. A. Khan, Y.-K. Park, K. K. Yadav, R. Vaziri, J.-H. Hwang, W. H. Lee, Y.-T. Ahn and B.-H. Jeon, Rev. Environ. Sci. Biotechnol., 2023, 22, 679–714.
8. D. B. Tripathy, ACS EST Water, 2024, 4, 4721–4740.
9. K. Gasljevic and E. Matthys, Energy Build., 1993, 20, 45–56.
10. H. Shi, W. Ge, H. Oh, S. M. Pattison, J. T. Huggins, Y. Talmon, D. J. Hart, S. R. Raghavan and J. L. Zakin, Langmuir, 2013, 29, 102–109.
11. M. Kotenko, H. Oskarsson, C. Bojesen and M. P. Nielsen, Energy, 2019, 178, 72–78.
12. G. Rose and K. Foster, J. Non-Newtonian Fluid Mech., 1989, 31, 59–85.
13. K. Gasljevic, G. Aguilar and E. Matthys, J. Non-Newtonian Fluid Mech., 2001, 96, 405–425.
14. J. Yang, Curr. Opin. Colloid Interface Sci., 2002, 7, 276–281.
15. R. K. Rodrigues, L. A. S. Silva, G. G. Vargas and B. V. Loureiro, J. Surfactants Deterg., 2020, 23, 21–40.
16. R. S. Mitishita, G. J. Elfring and I. A. Frigaard, J. Non-Newtonian Fluid Mech., 2022, 301, 104724.
17. R. R. Gupta, M. Daneshi, I. Frigaard and G. Elfring, Soft. Matter., 2024, 20, 4715–4733.
18. Z. Chu, C. A. Dreiss and Y. Feng, Chem. Soc. Rev., 2013, 42, 7174–7203.
19. H. Yu, Z. Xu, D. Wang, X. Chen, Z. Zhang, Q. Yin and Y. Li, Polym. Chem., 2013, 4, 5052.
20. C. Lim, J. D. Ramsey, D. Hwang, S. C. M. Teixeira, C. Poon, J. D. Strauss, E. P. Rosen, M. Sokolsky-Papkov and A. V. Kabanov, Small, 2022, 18, 2103552.
21. W. Hartt, L. Bacca and E. Tozzi, Soft Matter Series, Royal Society of Chemistry, Cambridge, 2017, pp. 379–398.
22. F. Choi, R. Chen and E. J. Acosta, J. Colloid Interface Sci., 2020, 564, 216–229.
23. J.-B. Salmon, A. Colin, S. Manneville and F. Molino, Phys. Rev. Lett., 2003, 90, 228303.
24. S. Lerouge and J.-F. Berret, Polymer Characterization, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, vol. 230, pp. 1–71.
25. S. M. Fielding and H. J. Wilson, J. Non-Newtonian Fluid Mech., 2010, 165, 196–202.
26. M. A. Fardin and S. Lerouge, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 91.
27. T. Divoux, M. A. Fardin, S. Manneville and S. Lerouge, Annu. Rev. Fluid Mech., 2016, 48, 81–103.
28. P. Rassolov, A. Scigliani and H. Mohammadigoushki, Soft. Matter., 2017, 13, 8368–8378.
29. L. Bécu, D. Anache, S. Manneville and A. Colin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 011503.
30. M. A. Fardin, D. Lopez, J. Croso, G. Grégoire, O. Cardoso, G. H. McKinley and S. Lerouge, Phys. Rev. Lett., 2010, 104, 178303.
31. M. A. Fardin, T. J. Ober, C. Gay, G. Grégoire, G. H. McKinley and S. Lerouge, EPL, 2011, 96, 44004.
32. J. Beaumont, N. Louvet, T. Divoux, M.-A. Fardin, H. Bodiguel, S. Lerouge, S. Manneville and A. Colin, Soft. Matter., 2013, 9, 735–749.
33. S. M. Fielding, J. Rheol., 2016, 60, 821–834.
34. R. G. Larson, E. S. G. Shaqfeh and S. J. Muller, J. Fluid Mech., 1990, 218, 573–600.
35. R. Larson, S. Muller and E. Shaqfeh, J. Non-Newtonian Fluid Mech., 1994, 51, 195–225.
36. P. Pakdel and G. H. McKinley, Phys. Rev. Lett., 1996, 77, 2459–2462.
37. G. H. McKinley, P. Pakdel and A. Öztekin, J. Non-Newtonian Fluid Mech., 1996, 67, 19–47.
38. 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, Phys. Rev. Fluids, 2022, 7, 080701.
39. C. A. Browne and S. S. Datta, Sci. Adv., 2021, 7, eabj2619.
40. C. A. Browne and S. S. Datta, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2320962121.
41. P. F. Salipante, C. A. E. Little and S. D. Hudson, Phys. Rev. Fluids, 2017, 2, 033302.
42. S. J. Haward, F. J. Galindo-Rosales, P. Ballesta and M. A. Alves, Appl. Phys. Lett., 2014, 104, 124101.
43. M. Y. Hwang, H. Mohammadigoushki and S. J. Muller, Phys. Rev. Fluids, 2017, 2, 043303.
44. Y. Zhang, H. Mohammadigoushki, M. Y. Hwang and S. J. Muller, Phys. Rev. Fluids, 2018, 3, 093301.
45. P. F. Salipante, S. E. Meek and S. D. Hudson, Soft. Matter., 2018, 14, 9020–9035.
46. R. M. Matos, M. A. Alves and F. T. Pinho, Exp. Fluids, 2019, 60, 145.
47. F. Hillebrand, S. Varchanis, C. C. Hopkins, S. J. Haward and A. Q. Shen, Soft. Matter., 2024, 20, 7133–7146.
48. G. R. Moss and J. P. Rothstein, J. Non-Newtonian Fluid Mech., 2010, 165, 1505–1515.
49. S. J. Haward, N. Kitajima, K. Toda-Peters, T. Takahashi and A. Q. Shen, Soft. Matter., 2019, 15, 1927–1941.
50. C. C. Hopkins, S. J. Haward and A. Q. Shen, Small, 2020, 16, 1903872.
51. C. Hopkins, A. Shen and S. Haward, Soft. Matter., 2022, 18, 8856–8866.
52. C. Hopkins, S. Haward and A. Shen, Soft. Matter., 2022, 18, 4868–4880.
53. S. Chen and J. P. Rothstein, J. Non-Newtonian Fluid Mech., 2004, 116, 205–234.
54. S. Wu and H. Mohammadigoushki, J. Rheol., 2018, 62, 1061–1069.
55. S. Wu and H. Mohammadigoushki, Soft. Matter., 2021, 17, 4395–4406.
56. G. R. Moss and J. P. Rothstein, J. Non-Newtonian Fluid Mech., 2010, 165, 1–13.
57. S. De, S. P. Koesen, R. V. Maitri, M. Golombok, J. T. Padding and J. F. M. van Santvoort, AIChE J., 2018, 64, 773–781.
58. S. J. Haward, C. C. Hopkins and A. Q. Shen, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2111651118.
59. A. Groisman and V. Steinberg, Nature, 2001, 410, 905–908.
60. T. Burghelea, E. Segre, I. Bar-Joseph, A. Groisman and V. Steinberg, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 066305.
61. J. Zilz, R. J. Poole, M. A. Alves, D. Bartolo, B. Levaché and A. Lindner, J. Fluid Mech., 2012, 712, 203–218.
62. L. Casanellas, M. Alves, R. Poole, S. Lerouge and A. Lindner, Soft. Matter., 2016, 12, 6167–6175.
63. A. Souliès, J. Aubril, C. Castelain and T. Burghelea, Phys. Fluids, 2017, 29, 083102.
64. P. Shakeri, M. Jung and R. Seemann, Phys. Fluids, 2022, 34, 073112.
65. H. Rehage and H. Hoffmann, J. Phys. Chem., 1988, 92, 4712–4719.
66. R. W. Mair and P. T. Callaghan, Europhys. Lett., 1996, 36, 719–724.
67. C. Oelschlaeger, M. Schopferer, F. Scheffold and N. Willenbacher, Langmuir, 2009, 25, 716–723.
68. S. J. Haward and G. H. McKinley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 031502.
69. S. J. Haward, T. J. Ober, M. S. Oliveira, M. A. Alves and G. H. McKinley, Soft. Matter., 2012, 8, 536–555.
70. D. Gaudino, R. Pasquino and N. Grizzuti, J. Rheol., 2015, 59, 1363–1375.
71. C. W. Macosko, Rheology: principles, measurements, and applications, VCH, New York, 1994.
72. Y. Zhao, A. Shen and S. Haward, Soft. Matter., 2016, 12, 8666–8681.
73. R. G. Larson, J. Rheol., 2012, 56, 1363–1374.
74. R. Granek and M. E. Cates, J. Chem. Phys., 1992, 96, 4758–4767.
75. M. S. Turner and M. E. Cates, Langmuir, 1991, 7, 1590–1594.
76. Y. Zhao, P. Cheung and A. Q. Shen, Adv. Colloid Interface Sci., 2014, 211, 34–46.
77. J. P. Rothstein and H. Mohammadigoushki, J. Non-Newtonian Fluid Mech., 2020, 285, 104382.
78. M. D. Abràmoff, P. J. Magalhães and S. J. Ram, Biophotonics Intern., 2004, 11, 36–42.
79. J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak and A. Cardona, Nat. Methods, 2012, 9, 676–682.
80. W. Thielicke and E. J. Stamhuis, J. Open Res. Softw., 2014, 2, e30.
81. R. J. Adrian and J. Westerweel, Particle Image Velocimetry, Cambridge University Press, New York, 2011.
82. M. Raffel, C. E. Willert, F. Scarano, C. J. Kähler, S. T. Wereley and J. Kompenhans, Particle Image Velocimetry: A Practical Guide, Springer, Switzerland, 2018.
83. L. Ducloué, L. Casanellas, S. Haward, R. Poole, M. Alves, S. Lerouge, A. Shen and A. Lindner, Microfluid. Nanofluid., 2019, 23, 33.
84. R. Rusconi, S. Lecuyer, N. Autrusson, L. Guglielmini and H. Stone, Biophys. J., 2011, 100, 1392–1399.
85. R. Poole, A. Lindner and M. Alves, J. Non-Newtonian Fluid Mech., 2013, 201, 10–16.
86. M. K. Kim, K. Drescher, O. Shun Pak, B. L. Bassler and H. A. Stone, New J. Phys., 2014, 16, 065024.
87. A. Groisman and V. Steinberg, New J. Phys., 2004, 6, 29.
88. D. W. Carlson, K. Toda-Peters, A. Q. Shen and S. J. Haward, J. Fluid Mech., 2022, 950, A36.
89. E. Y. Chen, C. A. Browne, S. J. Haward, A. Q. Shen and S. S. Datta, Stagnation points at grain contacts generate an elastic flow instability in 3D porous media, arXiv, 2024, preprint, arXiv:2412.03510 [physics] DOI:10.48550/arXiv.2412.03510,https://arxiv.org/abs/2412.03510.
90. C. Browne, A. Shih and S. Datta, J. Fluid Mech., 2020, 890, A2.
91. M. Kumar, S. Aramideh, C. A. Browne, S. S. Datta and A. M. Ardekani, Phys. Rev. Fluids, 2021, 6, 033304.
92. E. Y. Chen and S. S. Datta, J. Rheol., 2025, 69, 235–249.
93. G. K. Batchelor, J. Fluid Mech., 1971, 46, 813–829.
94. D. V. Boger, Annu. Rev. Fluid Mech., 1987, 19, 157–182.
95. A. Mongruel and M. Cloitre, Phys. Fluids, 1995, 7, 2546–2552.
96. A. Mongruel and M. Cloitre, J. Non-Newtonian Fluid Mech., 2003, 110, 27–43.
97. L. Rodd, J. Cooper-White, D. Boger and G. McKinley, J. Non-Newtonian Fluid Mech., 2007, 143, 170–191.
98. P. Coussot, J. Non-Newtonian Fluid Mech., 2014, 211, 31–49.
99. N. Waisbord, N. Stoop, D. M. Walkama, J. Dunkel and J. S. Guasto, Phys. Rev. Fluids, 2019, 4, 063303.
100. M. S. Abdelgawad, S. J. Haward, A. Q. Shen and M. E. Rosti, Phys. Fluids, 2024, 36, 113127.
101. D. Kawale, G. Bouwman, S. Sachdev, P. L. J. Zitha, M. T. Kreutzer, W. R. Rossen and P. E. Boukany, Soft. Matter., 2017, 13, 8745–8755.
102. Y. Ichikawa and M. Motosuke, Microfluid. Nanofluid., 2022, 26, 44.
103. M. Chellamuthu and J. P. Rothstein, J. Rheol., 2008, 52, 865–884.
104. S. A. Rogers, M. A. Calabrese and N. J. Wagner, Curr. Opin. Colloid Interface Sci., 2014, 19, 530–535.

---

Footnote

† Electronic supplementary information (ESI) available: Movie S1. Steady laminar base flow in the R_i = 300 μm channel at Wi = 0.18 < Wi_c. Successive half-loops are imaged with a time offset of 2 min. Scale bar is 500 μm. Movies play at 3.3× real time and are generated using a moving average of fluorescence intensity over 30 frames. Movie S2. Strongly 3D unsteady flow in the R_i = 300 μm channel at Wi = 0.78 > Wi_c. Successive half-loops have a time offset of 2 min. Scale bar is 500 μm. Movies play at 3.3× real time and are generated using a moving average of fluorescence intensity over 30 frames. Movie S3. Dead zones arise in the unsteady flow regime in serpentine channels of varied radii of curvature: R_i = 300, 500, and 1000 μm channels, from left to right. Flow is from left to right. Scale bar is 500 μm. Movies play at 3.3× real time and are generated using a moving average of fluorescence intensity over 30 frames. Movie S4. 3D twisting events arise in the unsteady flow regime in serpentine channels of varied radii of curvature: R_i = 300, 500, and 1000 μm channels, from left to right. Flow is from left to right. Scale bar is 500 μm. Movies play at 3.3× real time and are generated using a moving average of fluorescence intensity over 30 frames. See DOI: https://doi.org/10.1039/d5sm00344j

---