Transient flow behavior of complex fluids in microfluidic channels

The ongoing development of microfluidic devices involves the use of highly complex fluids, even of multiphase systems. Despite the great achievements in the development of numerous applications, there is still a lack in the complete understanding of the underlying physics of the observed macroscopic effects. One prominent example is the flow through benchmark contractions where micro- and even macroscopic explanations of some of the occurring flow patterns are still missing. Here, we study the development of the flow profiles of shear thinning semi-dilute polymer solutions in microfluidic planar abrupt contraction geometries. Flow profiles along the narrow channel part are obtained by μ-PIV measurements, whereby the pressure drop along the microfluidic channel as well as the local transient viscosities downstream to the orifice are computed. A relaxation process of the flow profiles from an initially parabolic shape to the flattened steady-state flow profile is observed and traced back to the polymer relaxation.


Introduction
Flow of complex fluids is ubiquitous to many applications of microfluidics (Stone et al. 2004). The resulting flow profiles depend on the nonlinear behavior of the medium and the channel geometry (van Steijn et al. 2007;Teh et al. 2008). While for continuous flow of complex fluids approximate solutions can predict the resulting flow profiles for standard geometries (Nguyen and Wereley 2002;Bruus 2009), much less is known on the phenomena occurring at sudden changes of channel geometry. One prominent example is the benchmark problem of entry flows into sudden contractions (Wassner et al. 1999;Alves et al. 2003;Lanzaro and Yuan 2011;Li et al. 2011;Hassell et al. 2011;Tamaddon Jahromi and Webster 2011;Sousa et al. 2011). Investigations on the entrance length of the flow of Newtonian fluids in microchannels (Ahmad and Hassan 2010) show that a steady-state flow profile is developed at a finite distance downstream to an orifice after the occurrence of small turbulences directly at the inlet. Nevertheless, with regards to non-Newtonian fluids, this situation is further complicated by the constriction geometry (Haward et al. 2010) as the extensional viscosity of the polymer dominates the fluid's response upon channel entry. Thus, the history of polymer flow prior to channel entrance has to be taken into account (Pérez-Trejo et al. 1996). Recent studies show that the length of a constriction has even impact on the upstream flow kinematics (Rodd et al. 2010). Although detailed studies on a single-molecule level have proven to be powerful in addressing the effect of elongation on the polymer configuration (Perkins 1997;Smith 1998;Larson et al. 1999;Squires and Quake 2005;Steinhauser et al. 2012) and efforts are made in simulating and calculating the flows of complex fluids on different length scales (Feigl 1996;Doyle et al. 1998;Xue et al. 1998;Lindner et al. 2003;Kröger 2004;Lubansky et al. 2007;Castillo-Tejas et al. 2009), microscopic explanations for the observable macroscopic flow effects remain challenging. To this end, detailed studies on the flow behavior of these fluids in benchmark geometries and its relation to microscopic mechanisms are essential.
Here, we study the flow of semi-dilute polymer solutions with non-Newtonian properties in microfluidic channels. We observe that the relaxation process of the flow profiles after a planar abrupt contraction takes place on a millimeter length scale and on a 0.5-5 s timescale. The flow profile relaxation time collapses to a single curve-by normalization to the characteristic relaxation time of the fluid-and decays with increasing Weissenberg number. We monitor the relaxation process by observing changes in the flow profile of the solution from an initially parabolic shape at the entrance to the expected flattened steady-state profile further downstream. The evolution of the flow profile is reproduced for all measured solutions.

Materials and methods
Microfluidic PDMS (Sylgard 184, World Precision Instruments) chambers with different geometries were built by the standard protocol (Duffy et al. 1998). A sketch of a geometry used for the measurements is shown in Fig. 1. For a fixed channel height of 80 lm, two constriction ratios of 7:1 (400:60 lm) and 10:1 (600:60 lm) were assembled. As the results of the two different constriction setups do not vary significantly, the reported measurements are not labeled to reflect the corresponding ratios.
The rheological properties of the solutions were determined by measurements using a rheometer with a cone-plate geometry (cone 1°0 0 25 00 , 40 mm, AR-G2, TA-Instruments). To account for possible aging effects of the fluids, samples of the polymer solutions (including the tracer particles) used for the microfluidic measurements were concurrently measured with the rheometer. Additional measurements of the PAA solutions without tracer particles showed that the particles did not affect the rheology outcome of the experimental solutions. Random measurements of the samples that had passed the microfluidic channel confirmed that the rheological properties of the polymer solutions did not change due to the channel flow. The measurement data obtained from the rheometer are best described by an adaptation of the Carreau model: c denotes the applied shear rate, g the shear viscosity, g 0 the zeroshear viscosity, and m eta a fit parameter describing the high shear rate behavior of the viscosity. The critical shear rate _ c c indicates the transition from the Newtonian (shear rate independent) viscosity plateau to the shear thinning region. From the critical shear rate we compute the characteristic relaxation time s rheo ¼ 1=_ c c of the fluid.
To vary the characteristic relaxation times of the sample fluids, the solvent viscosity was changed by replacing water with glycerol in varying volume ratios (ranging from 0 to 50 %). The final PAA concentration was adjusted to 2 % (w/v). Unless otherwise denoted, the various sample fluids are marked with the same colors throughout the figures, and are noted in Table 1.
The measurement setup consists of a Zeiss Axiovert 100 TV microscope with a 639 long distance objective (NA 0.75, resulting in a depth of field of 6 lm) and a motorized microscopy stage (Thorlabs Max 203) on which the microfluidic channel was mounted.
A syringe pump (SP210IWZ, World Precision Instruments) was used to introduce a volume-controlled flow. The pressure at the inlet of the microfluidic channel was continually measured with a pressure transducer (Sensortechnics, 24PC02K0xxA) to ensure constant flow rates.
Videos were captured at different positions along the microfluidic channel using either a Hamamatsu Orca or an  Andor Neo SCMOS camera with OpenBox (Schilling et al. 2004) or Andor Solis Software, respectively, for recording. The field of view was on the order of 100 lm 9 100 lm and the measurement plane was chosen to be at the vertical center of the channel (Fig. 1). To cancel out any timedependent measurement effects (e.g., heating of the microfluidic chamber, drift), the recording positions of the video-samples along the flow direction were not unidirectional or contiguous. The first measurement position at x = 0 mm was chosen such that the slightly round edges of the channel constriction (with a radius of *5 lm) were just out of view. The videos were analyzed by in-house developed Matlab scripts. The application of a brightness threshold filter to the video images reduced the depth of field to about 3 lm. To improve the measurement accuracy, a software filter for particle aggregates was also applied. A standard PIV algorithm (Williams et al. 2010) was implemented and applied to determine the flow profiles (v x y ð Þ). 1,000-3,000 image pairs were correlated for each flow profile to achieve high accuracy.
The obtained flow profiles are fitted by a power law v x y ð Þ ¼ ajy À bj ð Þ mþ1 þc. Thereby a is a parameter taking into account the volume flux, c the maximum flow speed and b the experimental center of the flow profile. The measure m characterizes how much the flow profile deviates from a parabolic flow profile (where m = 1). This fit is used as it is an exact solution for the flow of power law fluids (g PL _ c ð Þ ¼ g 1 _ c 1=m eta À1 ) in a circular tube or infinite parallel plate geometry with m ¼ m eta (Bruus 2009). Due to the geometric boundary conditions the parameter m, obtained by fitting the fully developed flow profile of a power law fluid in rectangular geometries, varies from m eta as measured by a rheometer. The correction factor is determined by numerical calculations, which show that in a rectangular geometry with the chosen aspect ratio of width to height of 0.75 we obtain m ¼ À0:24 þ 1:32 m eta .
To eliminate the artifacts due to beads sticking to the micro channel walls and roughness of the channel walls, only data points with v [ eÁc and correlation values above a certain threshold were used for the power law profile fit. To evaluate the pressure drop and the transient viscosities (Figs. 3 and 4) e was set to 0.2; for all other measurements, we set e = 0.7 to simplify the fitting procedure. The change in the value of e leads to only slight changes in the obtained fit parameters (±10 %).
To evaluate the relaxation of the fluid after a constriction (Fig. 1), a series of measurements of a single polymer solution at a constant flow rate was done at different positions. The resulting temporal evolution of the flow profiles was fit with the equation where m 0 denotes m at x = 0, m end the steady-state value, and l ss the decay length. The single exponential decay function was the simplest model, which fit the data bestindicating that there is one dominating relaxation process governing the response of the fluid after the contraction.
To check the flow conditions in the wider part of the channel, velocity profiles in both directions perpendicular to the flow direction were recorded far upstream of the contraction such that effects of the inlet and the orifice were negligible. The shape of the profile along the shorter axis showed the expected power law profile with the characteristic flattening, whereas the profile along the longer axis showed a pronounced ''rounding''. These Data presented in Figs. 2b, 3, and 4 are obtained from the same measurement series of fluid F3 (Table 1) at Wi ¼ 290.

Results and discussion
The semi-dilute polymer solutions showed a pronounced shear thinning behavior, as determined by rheometer measurements (Colby et al. 2006) (Fig. 2a). Their respective characteristic relaxation times were determined to range from 0.9 to 23 s (Table 1). Flow profiles of such non-Newtonian fluids in rectangular flow chambers deviate from the flow profile expected for Newtonian fluids in similar three-dimensional geometries (Bruus 2009). While at the channel walls the shear rate is highest, the shear rate in the middle of the channel approaches zero. Thus, the shear thinning behavior of the viscosity occurs at the walls, while in the center of the channel higher viscosities dominate the flow. Consequently, the flow profile flattens in the middle of the channel (Fig. 2b).
All observed flow profiles are excellently described by a power law fit (Fig. 2b). In the limit of purely Newtonian fluids (g _ c ð Þ ¼ const), a parabolic flow profile is observed corresponding to m % m eta ¼ 1. For the non-Newtonian fluid F3 (Table 1), the parameter determined by a Carreau fit of m eta ¼ 2:1 should result in a power law profile in the microfluidic channel of m = 2.6. Indeed, in good agreement an exponent of m = 2.8 ± 0.1 was observed in steady-state conditions in the microfluidic channels (Fig. 3).  ). The inset shows the development of the local viscosity at different y-positions (y: green 1 lm, blue 2 lm, cyan 5 lm, red 10 lm, and black 25 lm). Near the center of the flow profile (y * 0), the viscosity increases with increasing distance from the orifice as the polymers, initially elongated due to the orifice flow, relax to their steady state in shear flow Once the fluid is forced into an orifice flow, the flow profile changes significantly. To quantify the development of the flow after the orifice, the fit exponent m at different measurement positions x along the narrow channel is plotted (Fig. 3). The fit of the velocity profile directly at the entrance of the narrow channel part at x = 0 mm yields m = 1. This corresponds to a parabolic flow profile. For 0 mm \ x \ 2 mm, the fit exponent shows a steep growth followed by a flattening until x = 4 mm, where it reaches a value of m = 2.8. The growth in the exponent m stems from a flattening in the center and simultaneous broadening of the flow profiles at the respective positions. This change in the flow profile is the dominating effect due to polymer relaxation. The volume constraint along the flow direction of the fluid causes a variation of the maximum flow speed of at most 15 %, which is also observed in the measurements. The steady-state value of m extracted from the flow measurements corresponds to the value found using coneplate shear rheometry, and is unchanged from x = 4 mm until the last measurement position at x = 20 mm. This relaxation behavior is replicated for all fluids studied.
At the orifice-the first measurement position-a parabolic shape is found for all Wi ¼ 10. . .1; 000 (Fig. 3 inset), indicating a constant local viscosity along the axis perpendicular to the flow direction. The effect of shear thinning on the steady-state flow profiles increases with increasing Wi from m = 2 to 4 (Fig. 3 inset). The reason is that the steady-state flow profile contains all shear rates from zero (at the center of the channel) to the respective shear rates at the walls. As the viscosity measurement shows a Newtonian plateau at low shear rates, the flow profile should still fit a parabola (m = 1) within a region near the center of the channel in contrast to the assumption of a power law velocity profile with m [ 1. With increasing flow speed, the shear rates increase and therefore the region in which the profile follows a parabola diminishes. This enhances the quality of the power law fit, since the fit is less influenced by the Newtonian viscosity plateau-the determinant of the flow profile center in steady state. Therefore, m end increases with increasing Wi.
The analysis of the shape of the flow profile allows for the calculation of the local pressure drop and hence an estimation of the local and transient viscosities at each measurement position. The calculation is based on the assumption that the contributions of elastic effects to the stress tensor are small compared to the viscous effects. This is justified, as the velocity components perpendicular (v y and v z ) to the flow direction (v x ) are negligible and are not observable by the PIV. Thus, the pressure drop is constant along the axes perpendicular to the flow direction at each measurement position x. The pressure drop can be calculated by assuming that the fluid elements near the walls are already in their steady state. This assumption can be justified by the flow conditions at the walls: the effective shear rate is always chosen such that Wi [ 1, thus, the shear flow dominates the polymer relaxation. Polymers near the wall experience the highest shear rate and thus reach their steady-state configuration faster (Takahashi et al. 1987). In addition, the flow speed near the channel wall is the lowest. Thus, at a given distance from the orifice polymers near the walls have enough time to relax to their steady-state configuration. This can be directly seen by computing the ratio of the corresponding time scales where t flow denotes the time a fluid element has travelled in the flow, b = y/w the fractional distance from the center of the channel, _ c w the corresponding shear rate and d the width of the recording area of a single video. For the given geometry layout and near the channel walls, we compute k ) 1 for all measurement positions. Thus, the polymers at the walls have travelled much longer in the fluid than the characteristic timescale of the flow (1=_ c w ). This further supports the notion that the polymers at the walls have reached their steady-state configuration by the first measurement position at x = 0 mm. From the law of conservation of momentum, it follows that with p denoting the pressure. In Fig. 4a, the pressure drop is plotted at different positions along the narrow channel. The pressure drop is lowest near the orifice. It increases from 1 Â 10 6 Pa m at the orifice to about 5 Â 10 6 Pa m with a decay length of ca. 1 mm. The steady state is reached at a distance of 4 mm from the orifice.
The calculated pressure drop is in good accordance with the overall pressure drop measured by a pressure transducer. This is expected, as the narrow microfluidic channel contributes most to the overall pressure drop between the pressure transducer and the outlet of the microfluidic channel at ambient pressure. Now the local transient viscosities g þ perpendicular to the flow direction (Fig. 4b) can be computed at each measurement position from the pressure drop using To prevent the calculation of infeasible outcomes near the center of the channel, where the shear rates fall to zero, the transient viscosities are limited to the corresponding steady-state viscosity values determined by the rheometer measurement.
It can be seen that directly at the orifice (x *10 lm) the viscosity is dramatically lower than for the measurement positions further downstream. Moreover, as this profile shows a nearly parabolic shape (m *1), the local viscosity perpendicular to the flow direction is almost constant. With increasing distance from the orifice, the shear thinning behavior, which manifests in the sharp dependence of the viscosity on the position perpendicular to the flow direction, becomes more pronounced. Near the middle of the channel-at y = 1 lm-a sharp increase in the viscosity is observed, with maximal viscosity attained at x * 4 mm (Fig. 4b inset). From there on, the local viscosity profile along the flow direction remains unchanged and the steady state is reached. Near the channel walls, the local viscosity shows a slight decrease. As the shear rates near the center of the channel decrease, observed by a flattening of the flow profiles, the shear rates near the channel walls increase to maintain the volume flux. Since the local viscosity at the channel walls is always in its steady state described by the Carreau model, the viscosity decreases moderately with an increasing distance from the orifice.
To better quantify the development of the flow profiles, the parameter m(x) is fitted with an exponential function (Fig. 3) and the obtained decay lengths l ss are plotted in Fig. 5 (inset) as a function of the dimensionless Weissenberg number Wi. The colors correspond to different test fluids with varying relaxation times ( Table 1). The decay length increases up to a factor of 7 for all measured solutions with increasing Wi (Fig. 5 inset). Already for moderate Wi of about 1-10, the decay length l ss is in the range of 0.3-1 mm; for high Wi ¼ 100. . .1; 000, l ss reaches up to 2 mm. Thus, the relaxation process plays an important role on a length scale that is much larger than the smallest channel dimension (10-100 lm). This is contrary to the common steady-state conditions in microfluidic devices where only the smallest dimensions of the channel must be considered.
To relate the observed relaxation process to the physical mechanism, we compute the decay time s channel from the decay length l ss by considering the maximal flow velocity at y = 0. Fig. 5 shows the normalized relaxation times of the flow profiles s norm ¼ s channel =s rheo in the microfluidic channel at various dimensionless Weissenberg numbers Wi, at varying flow speeds and for fluids with different characteristic relaxation times. This results in a collapse of s norm to a single curve, which decays with increasing Wi. For moderate Wi $ 10, the measured relaxation time in the microfluidic channel equals the characteristic relaxation time of the polymer solutions (s norm ¼ 1). With increasing flow speed, the relaxation time of the flow profile in the microfluidic channel is drastically decreased to about 0:04s rheo at Wi ¼ 1; 000.

Conclusion
While the flow of Newtonian fluids in microfluidic setups always obeys parabolic flow profiles, the non-Newtonian behavior of polymer solutions introduces a history dependence, which needs to be taken into account. Shear thinning fluids in channel flows show a flattened flow profile, which is best described by a power law, whereby varying flow geometries complicate matters further. Here, we present the effect of an extensional flow geometry, where we observed a parabolic flow profile of a semi-dilute polymer solution directly after an orifice, and that appears to be independent of the flow rate. The subsequent relaxation process is flow rate dependent, yet can be renormalized by the characteristic relaxation time of the fluid.
The parabolic shape of the flow profiles measured directly after the orifice can be explained by considering the flow conditions directly upstream of the orifice. The extensional flow caused by the orifice leads to an extension of the polymers. Consequently, the resistance of the polymers to the extensional flow increases and results in an adjustment of the extensional flow to minimize the dissipated energy. Finally, all polymers along the axis perpendicular to the flow direction are equally extended. As all polymers along this axis have similar average configurations, they affect the same resistance to the flow. Thus, a parabolic flow profile is observed. Above a certain extensional rate, which is mandatory to start the extension of the polymers from their equilibrium conformation, this outcome is independent of the flow rate (Perkins 1997;Smith 1998;Larson et al. 1999). The relaxation of the initially stretched polymers in shear flow after the contraction causes an increase in their hydrodynamic diameter and consequently an increase in the dissipation of energy  (Colby et al. 2006). Therefore, the local viscosities along the flow direction increase until the polymers reach their steady-state conformation.
For moderate Wi $ 10, the time that the fluid needs to reach its steady-state flow profile equals the characteristic relaxation time. At such low Wi, the shear rate in the center of the channel (-6 lm \ y \ 6 lm) is lower than the critical shear rate. Therefore, the polymers relax to their equilibrium conformation within their characteristic relaxation time s rheo .
With increasing Wi, the relaxation time of the flow profile decreases. The reason is that at high Wi, the region around the center of the flow profile, where the local shear rate is smaller than 1=s rheo , diminishes and thus only marginally affects the overall flow profile. Consequently, the flow profile is determined by the relaxation of the polymers at high shear rates _ c [ 1=s rheo . In this regime, polymers that are exposed to higher shear rates reach their steady state faster, since with higher shear rates only modes of the polymers with short relaxation times can relax (Colby et al. 2006). Consequently, the steady state of the polymer configuration is further away from the initially stretched state for polymers exposed to low shear rates than for polymers exposed to higher shear rates. Therefore, this relaxation is observable for the longest time near the symmetry axes of the channel (y * 0) where it causes a flattening of the flow profile and the relaxation time of the complete flow profile decreases with increasing Wi.
The detailed insights of the flow profile relaxation of polymer solutions may pave the way for a broader development of complex microfluidic devices. Generally, the physical phenomena presented in this work play a role in all flows of complex fluids in microfluidic channels with extensional flow geometries. In particular, fluids in which anisotropies can be induced by such geometries may show the relaxation process examined here. Many of the commonly used fluids, such as emulsions and suspensions containing large biomolecules, industrial polymers or other non-spherical particles, will show the presented effects. The long relaxation times may turn out to be important for bio-analytical assays. There may also be ways of harnessing such local viscosity differences in chemical reactions in emerging lab-on-a-chip applications.