Flow-based manipulation of particles is an essential tool for studying soft materials, but prior work has nearly exclusively relied on using two-dimensional (2D) flows generated in planar microfluidic geometries. In this work, we demonstrate 3D trapping and manipulation of freely suspended particles, droplets, and giant unilamellar vesicles in 3D flow fields using automated flow control. Three-dimensional flow fields including uniaxial extension and biaxial extension are generated in 3D-printed fluidic devices combined with active feedback control for particle manipulation in 3D. Flow fields are characterized using particle tracking velocimetry complemented by finite-element simulations for all flow geometries. Single colloidal particles (3.4  μm diameter) are confined in low viscosity solvent (1.0 mPa s) near the stagnation points of uniaxial and biaxial extensional flow for long times ( 10 min) using active feedback control. Trap stiffness is experimentally determined by analyzing the power spectral density of particle position fluctuations. We further demonstrate precise manipulation of colloidal particles along user-defined trajectories in three dimensions using automated flow control. Newtonian liquid droplets and GUVs are trapped and deformed in precisely controlled uniaxial and biaxial extensional flows, which is a new demonstration for 3D flow fields. Overall, this work extends flow-based manipulation of particles and droplets to three dimensions, thereby enabling quantitative analysis of colloids and soft materials in complex nonequilibrium flows.

The ability to trap and manipulate single particles in liquid media has revolutionized several fields of science and engineering. Particle trapping methods based on optical [1], magnetic [2], acoustic [3], and electrokinetic [4] force fields have allowed for detailed investigations of colloids and soft materials. Hydrodynamic flow-based traps have further enabled the study of freely suspended particles such as liquid droplets, single polymers, and vesicles in nonequilibrium extensional flow conditions [5–8]. In recent years, the Stokes trap has allowed for multiplexed particle trapping and manipulation using automated flow control, including the fluidic-directed assembly of colloidal particles [9]. The Stokes trap relies on flow-based manipulation of freely suspended particles or molecules to study the conformation, deformation, and dynamics of materials using a model predictive control (MPC) algorithm [9]. Previous generations of the flow-based trapping method allowed for the confinement of single particles or molecules near the stagnation point of planar extensional flow generated in polydimethylsiloxane (PDMS)-based cross-slot microfluidic devices fabricated using standard soft lithography [10–12].

To date, flow-based trapping has nearly exclusively relied on using a planar microfludic device and model approximation to manipulate the 2D center-of-mass positions of particles [7–9,13]. Prior generations of flow-based traps relied on 2D manipulation in thin-gap (or Hele–Shaw) microdevices due to the limitations imposed by the layer-by-layer fabrication method used in soft lithography [9–12,14]. However, understanding the dynamics of soft materials (e.g., polymeric and inkjet solutions and colloidal clusters) in 3D flows is highly relevant for industrial processes that rely on uniaxial extensional and biaxial extensional flow [15–18]. For example, electrospinning and extrusion processes generally subject polymeric liquids to uniaxial extensional flow fields [19–21], whereas film blowing and vacuum molding processes subject materials to biaxial extension [22,23]. Moreover, 2D flows such as planar extension have been reported to give rise to different apparent extensional viscosities for viscoelastic liquids compared to 3D uniaxial and biaxial extensional flows [24]. Broadly speaking, polymeric liquids are commonly subjected to processing and deformation in 3D extensional flows in a wide range of industrial processes [25], but generating these 3D flows in micro- or millifluidic devices is challenging and relatively uncommon [26]. Prior work has reported extensional viscosity measurements in uniaxial and biaxial extensional flows using the opposed-jets device, typically with millimeter nozzle dimensions or larger [27–29].

Microfluidic devices capable of generating 3D flows have been reported based on a 3D stacked, multilayer six-arm cross-slot geometry [30]. Here, Gonzalez and Liu used computational modeling to characterize fluid flow in a double-layered microfluidic design with the overall goal of translating particles in 3D using cooperative flow modes along the three spatial axes [30]. By operating this device in an opposing flow mode wherein channels on opposite sides serve as inlets, a flow field with 3D hyperbolic streamlines was generated accompanied by a stagnation point near the center of the device [30]. Prior work has also studied elastic instabilities in six-arm [31] and four-arm [32,33] channel geometries with finite aspect ratios using computational modeling. However, creating highly structured 3D geometries in PDMS-based devices is challenging due to device fabrication constraints associated with channel dimensions and flow channel aspect ratios [14]. Haward et al. recently reported a six-arm cross-slot geometry device fabricated by selective laser-induced etching in quartz substrates [26,34]. Uniaxial and biaxial extensional flows were generated in this device, and the birefringence signals of dilute polymer solutions were measured as a function of flow strength [26,34].

In recent years, additive manufacturing has provided a powerful alternative to PDMS-based soft lithography techniques to fabricate 3D microfluidic devices using a variety of materials [35]. Additive manufacturing allows for fabrication of well-defined 3D flow channels, enabling researchers to build complex shapes as well as nonorthogonal and nonplanar structures inside microfluidic systems. Using stereolithography (SLA), Lee et al. fabricated a 3D printed microfluidic device with a helical microchannel to detect bacterial cells [36]. In addition, Spivey et al. created a single-cell capturing device using projection printing (DMD-PP) to observe cellular aging [37]. Despite recent progress, however, automated flow control has not been extended to trap and manipulate particles in 3D flows.

In this work, we demonstrate 3D trapping and manipulation of freely suspended particles, droplets, and vesicles using automated flow control. SLA 3D-printing is used to fabricate flow devices with a six-arm cross-slot geometry [Figs. 1(a) and 1(b)]. Uniaxial and biaxial extensional flow fields are generated within 3D-printed devices using pressure-driven flow and are characterized using particle tracking velocimetry (PTV) and finite-element simulations. Using automated flow control, single colloidal particles are trapped and confined near the stagnation point in both uniaxial and biaxial extensional flows, and trap performance is characterized by determining trap stiffness from Fourier spectral analysis of particle position fluctuations. We further demonstrate full 3D control over the center-of-mass position of trapped particles along user-defined trajectories using automated flow control. In addition, liquid droplets and giant unilamellar vesicles (GUVs) are deformed in well-defined uniaxial and biaxial extensional flows. Overall, this work significantly advances the ability to confine, manipulate, and study freely suspended particles in 3D flows, which will enable quantitative comparisons between microscopic dynamics to macroscopic rheological measurements.

FIG. 1.

Schematic of the 3D printed flow device used in this work. (a) Uniaxial and (b) biaxial extensional flow configurations using the six-channel cross-slot geometry. Red arrows denote inlets, whereas yellow arrows denote outlets. (c) Schematic (not to scale) of the fluidic device when mounted on an inverted microscope setup during autofocusing using the single-imaging mode scheme. Here, the microscope objective is fixed, and the device is translated orthogonal to the image plane (in the z-direction) to maintain particles in focus. Blue regions denote internal flow channels. The bottom channel is sealed with a glass coverslip. For clarity, the y-axis channels orthogonal to the plane of the page are not shown. (d) Schematic (not to scale) of orthogonal imaging setup wherein two objectives actively image and track a trapped particle, allowing for full 3D particle tracking. (e) Photograph of the 3D-printed device used in the single-imaging mode with internal channels filled with red dye. (f) Photograph of 3D-printed device used in dual-imaging mode with internal channels filled with blue dye. (g) Schematic of the overall device showing a cylinder supplying nitrogen gas to computer-controlled pressure regulators which in turn pressurize solution reservoirs. The reservoirs feed the solution into the single-imaging 3D-printed device (photograph shown in top-down view). Figure is shown in color in online version.

FIG. 1.

Schematic of the 3D printed flow device used in this work. (a) Uniaxial and (b) biaxial extensional flow configurations using the six-channel cross-slot geometry. Red arrows denote inlets, whereas yellow arrows denote outlets. (c) Schematic (not to scale) of the fluidic device when mounted on an inverted microscope setup during autofocusing using the single-imaging mode scheme. Here, the microscope objective is fixed, and the device is translated orthogonal to the image plane (in the z-direction) to maintain particles in focus. Blue regions denote internal flow channels. The bottom channel is sealed with a glass coverslip. For clarity, the y-axis channels orthogonal to the plane of the page are not shown. (d) Schematic (not to scale) of orthogonal imaging setup wherein two objectives actively image and track a trapped particle, allowing for full 3D particle tracking. (e) Photograph of the 3D-printed device used in the single-imaging mode with internal channels filled with red dye. (f) Photograph of 3D-printed device used in dual-imaging mode with internal channels filled with blue dye. (g) Schematic of the overall device showing a cylinder supplying nitrogen gas to computer-controlled pressure regulators which in turn pressurize solution reservoirs. The reservoirs feed the solution into the single-imaging 3D-printed device (photograph shown in top-down view). Figure is shown in color in online version.

Close modal

Devices are fabricated using an SLA 3D printer (Form 3, Formlabs) utilizing commercially available photopolymer resin (Clear V4, Formlabs). A 3D model of the device containing 750  μm-wide channels was sliced using the software Preform (Formlabs) and subsequently printed at a fixed 25  μm layer height. Devices were removed from print supports, and each internal channel within the device was copiously flushed with isopropyl alcohol via a syringe. Excess unpolymerized resin was washed away by submerging the devices in an isopropyl alcohol bath under strong agitation. The devices were then thoroughly dried, heated to 60  °C, and exposed to 405 nm light for 15 min to fully cure the 3D printed part according to manufacturer specifications. In some cases, mechanical sanding was necessary to smooth surfaces that retained print support features.

For single-mode imaging experiments involving one camera [Figs. 1(c) and 1(e)], the device design includes a channel with one open surface in the - z direction arm of the six-channel cross-slot. A glass coverslip was used in conjunction with UV-curing adhesive (Norland) to seal the bottom channel, thereby providing optical access to the center of the six-arm cross-slot. The distance from the coverslip on the bottom of the device to the center of the 3D cross-slot is 2 mm.

For orthogonal-imaging experiments involving two cameras [Figs. 1(d) and 1(f)], we utilized a different 3D-printed device with a modified design (vide infra). A 3D cross-slot geometry was positioned near a corner of the overall 3D-printed block such that two open surfaces on the + x and - z-direction channels were exposed. Glass coverslips encased these open channels using UV-curing adhesive as previously described, enabling imaging into the center of the 3D cross-slot from both the x y and y z planes. The distance from the two coverslip surfaces to the center of the cross-slot is 2 mm.

Numerical simulations were performed using COMSOL Multiphysics to model the flow of a Newtonian fluid in the 3D device geometry. A model of the printed 3D cross-slot was directly imported into COMSOL Multiphysics. For uniaxial extensional flow, four inlet channels were oriented along the y- and z-axes, and two outlet channels were oriented along the x-axis. For biaxial extensional flow, two inlet channels were oriented along the z-axis, and four outlet channels oriented along the x- and y-axes. The channel size of the simulated geometry was set to 750  μm to match the channel dimensions in experiments. In all cases, a no-slip boundary condition is specified for all internal surfaces, and static pressure values for outlet port cross sections and inlet port cross sections are imposed for steady-flow such that the inlet pressures are higher than outlet pressures. The geometry was meshed using a physics-controlled sequence and an extra fine mesh size, resulting in 246 204 elements. Finally, the velocity profile within the channels was solved using the Stokes equations.

Imaging for the single-camera microscope mode was performed using an inverted microscope (Olympus IX81) coupled with a CMOS camera (Grasshopper 3, FLIR) or an electron-multiplying charged-coupled device (EMCCD) camera (Andor iXon Ultra 897), which captured images at 1024 × 768 and 512 × 512 pixel resolutions, respectively. A 0.25 numerical aperture (NA) 10 ×, a 0.40 NA 10 ×, or a 1.10 NA 60 × long working distance objective (Olympus) was used. In some cases, an additional 1.6 × magnification lens was employed, resulting in a pixel size of 0.345 or 0.216  μm for the 10 × objectives and 0.267 or 0.167  μm for the 60 × objective. The 0.40 NA 10 × objective was used for particle manipulation experiments due to the smaller depth-of-field, which generally yielded enhanced z-position localization. The 1.10 NA 60 × objective was used for observing GUVs due to the higher NA, which aids in accurate and detailed observation of vesicles. Epi-illumination of the cross-slot was achieved using a 532 or a 635 nm continuous wave (CW) laser (Spectra-Physics) with a dichroic beam splitter (ZT532 or ZT650, Chroma). The emitted light passed through a 542 or a 655 nm long-pass filter (ET542lp or ET655lp, Chroma) prior to the camera.

For two-plane orthogonal imaging experiments, the dual-camera microscope mode consisted of the components from the single-imaging mode with an additional optical train consisting of a second CMOS camera (Grasshopper 3, FLIR), a tube lens (ThorLabs), and 542 nm long-pass filter (ET542lp, Chroma). The 10 × 0.25 NA and the 10 × 0.40 NA objectives were simultaneously utilized to generate two orthogonal imaging planes for fine-scale center-of-mass particle localization in 3D. The 0.40 NA objective ( x y plane objective) was used on the inverted microscope, and the 0.25 NA objective ( y z plane objective) was used in the second optical train.

Aqueous solutions (viscosity η = 1.0 mPa s) seeded with trace amounts ( < 1.8 × 10 3 beads/ μl) of 1.0  μm (Fluospheres) or 3.4  μm diameter (Spherotech) fluorescent Nile Red or Nile Blue polystyrene beads were introduced into the device using microfluidic tubing connected to a reservoir for the channels in the cross-slot geometry. Standard microfluidic fittings were used to interface the microfluidic tubing with Luer-lock ports embedded within the device design [Fig. 1(g)]. The pressures in each of the six reservoirs were individually controlled from the laboratory computer via either electropneumatic pressure regulators (Proportion-Air) using data acquisition cards (CompactDAQ, National Instruments) or piezoelectric pressure regulators (OB1 MK3, Elveflow). The latter set of pressure regulators were used in 3D orthogonal manipulation and droplet deformation experiments.

Uniaxial and biaxial extensional flow fields were generated near the center and axes of the cross-slot by pressurizing the inlet and outlet channels appropriate for the desired flow fields. Near the sidewalls of the channels, however, flow adjacent to the channel walls is dominated by shear due to the no-slip boundary condition. In this work, the horizontal imaging axis was chosen as the axis of extension ( x-axis) for uniaxial extensional flow. Biaxial extension was generated using the out-of-plane axis ( z-axis) as the compressional axis. For flow field characterization experiments, only select inlet channels were seeded with fluorescent beads to minimize background signal from out-of-focus beads. For uniaxial extensional flow, the + y and - y inlet channels contained fluorescent beads and the + z and z inlet channels contained pure solvent. For biaxial extensional flows, fluorescent beads were added to the + z channel. PTV was performed using previously reported algorithms to characterize the flow fields generated and to determine strain rates as a function of applied pressure (supplementary material, Figs. 1 and 2 [61]) [38].

A LabVIEW program was used for automated closed-loop feedback control of the imposed flow fields using proportional-integral-derivative (PID) control, as previously reported (supplementary material, Fig. 3 [61]) [11]. In the single-imaging mode, the particle position in the x y plane is determined by localizing center-of-mass positions from camera images. The pressure regulators are modulated via the PID controller to control imposed flow rates in the six channels. Two PID controllers are used in parallel to manipulate particle positions along the x- and y-axes by perturbing pressure gradients Δ p x and Δ p y, respectively. The total feedback loop time constant for the image acquisition in this work was 25–50 ms (40–20 Hz), though the current setup enables particles to be imaged and trapped using feedback rates up to 62 Hz.

In the single-imaging mode [Fig. 1(c)], particles are tracked in the z-axis direction (orthogonal to the image plane) by mounting the flow device on a motorized stage (Zaber Technologies) and translating the stage along the z-axis (supplementary material, Fig. 4 [61]). The velocity of the translational motor was controlled via LabVIEW, and the motor position was used to determine the z-position of trapped particles. An autofocusing algorithm was employed to maintain the particle of interest in the focal plane, as previously reported (supplementary material, Fig. 4 [61]) [39]. Using this approach, the z stage continually moved to track particles in the z-direction, and the average intensity I of the pixels corresponding to the detected particle was recorded as a function of the z position. The resulting fluorescence intensity gradient d I / d z was used to determine the translational direction of the z motor such that I was maximized, which corresponds to an estimated z position where the particle was in focus. In this work, we used the particle image intensity as the image quality parameter as opposed to the trapped particle area used by Hsiao et al. [39] due to the lower-magnification objective used in this work.

For particle manipulation experiments in 3D using the orthogonal dual-imaging mode [Fig. 1(d)], particles are directly imaged in the y z plane using an additional camera and optical imaging path, enabling direct determination of the z positions of trapped particles (supplementary material, Fig. 5 [61]). The dual-imaging approach provides greatly improved accuracy for determining the z-position of the trapped particle. In the dual-imaging mode, particle x, y, and z positions are acquired simultaneously from two imaging planes ( x y and y z planes). Three PID controllers are simultaneously employed to manipulate particle positions along the x-, y-, and z- axes (supplementary material, Fig. 6 [61]).

A solution of polydisperse water-in-oil droplets was prepared by repeatedly inverting a mixture containing 5  μl of 1.9 mM sodium dodecyl sulfate (SDS) in aqueous solution and 500  μl of 0.5 wt. % sorbitan monooleate surfactant (Span 80) in mineral oil. The combination of SDS and Span 80 has previously been reported to produce ultralow oil-water interfacial tensions γ 10 6 10 5 N/m [40–42]. However, to enable efficient delivery of droplets prepared off-device, slightly lower surfactant concentrations were used yielding an interfacial tension γ 10 4 10 3 N/m (as determined from droplet deformation experiments). For droplet deformation experiments, 0.5 wt. % Span 80 in mineral oil was added to each of the six reservoirs. The stock droplet solution was diluted 10-fold in the + z channel. A broadband LED mounted above the device was used to illuminate the 3D cross-slot, and droplets were imaged using brightfield microscopy. Trapped droplets were deformed under both uniaxial and biaxial extensional flow. The center-of-mass of a trapped droplet was tracked, and automated PID control was used to position the droplet near the stagnation point in the x- and y-directions. Movies of droplet deformation (supplementary material, Movies 12 and 13 [61]) were recorded as the pressure differentials were incrementally increased to observe drop deformation as a function of flow strengths. The edges of the trapped droplet were tracked using a custom LabVIEW program to quantify droplet deformation as a function of applied pressure.

GUVs were prepared from a mixture of 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC, Avanti Polar Lipids) and 0.2 mol. % of headgroup-labeled fluorescent lipid 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine-N-(Cyanine 5) (DOPE-Cy5, Avanti Polar Lipids) using a previously described electroformation method [43–45]. To prepare the stock lipid solution for electroformation of GUVs, a mixture of chloroform and methanol (in a 5:2 v/v ratio) was used. A small aliquot (7  μl) of the lipid mixture was then deposited onto indium tin oxide (ITO) coated glass slides as a thin film and dried under vacuum for 3 h. Next, 100 mM sucrose solution ( η = 1.1 mPa s) was introduced into the electroformation chamber, which was constructed using two of the prepared ITO slides and a PDMS spacer. The electroformation chamber was incubated overnight at 60  °C while connected to a function generator (Agilent 33 220 A) that applied an alternating current (AC) electric field of 10 V mm−1 at 10 Hz. To observe the GUVs using fluorescence microscopy, a 635 nm CW laser was used to illuminate the Cy5-labeled vesicles. Images were captured using an EMCCD camera at a frame rate of 18 Hz. Similar to the droplet stretching experiments, GUVs were trapped and deformed under both uniaxial and biaxial extensional flow using PID control.

We began by characterizing the 3D flow fields generated in the 3D-printed fluidic devices. Uniaxial extensional flow consists of two axes of compression and one axis of extension. Uniaxial extensional flow, therefore, has nonzero velocity in the z-axis and is described by
v x = ε ˙ u x ,
(1a)
v y = 1 2 ε ˙ u y ,
(1b)
v z = 1 2 ε ˙ u z ,
(1c)
where the origin of the coordinate system is located at the stagnation point and ε ˙ u is the strain rate in uniaxial extension. For uniaxial extensional flow, the x and y axes are defined as the two axes coincident with the imaging plane, and the z axis is orthogonal to the x y plane.
For comparison, planar extensional flow is a 2D flow field consisting of one axis of compression and an orthogonal axis of extension, as described by
v x = ε ˙ p x ,
(2a)
v y = ε ˙ p y ,
(2b)
v z = 0 ,
(2c)
where ε ˙ p denotes the strain rate for planar extensional flow.

Uniaxial extensional flow is generated by pressurizing the four inlet channels corresponding to the compressional y- and z-axes at equal positive values ( p y = p + y = p z = p + z) relative to the pressure in the two outlet reservoirs ( p x = p + x) along the extensional x-axis. We began by characterizing the flow using the single-mode imaging setup [Fig. 1(c)]. Figure 2(a) shows particle trajectories in uniaxial extensional flow obtained using PTV for imaging in the x y plane at an applied pressure Δ p u = p + z p + x = 1.4 kPa. Figure 2(a) also shows characteristic streamlines obtained from numerical simulations of the identical flow device geometry. Individual particles follow 3D streamlines and enter the focal plane of imaging (but are not necessarily in focus) due to nonzero v z (supplementary material, Movie 4 [61]). As expected, particle velocities v x and v y in the x y image plane are independent of z-position near the stagnation point and away from channel walls. The flow field is hyperbolic near the stagnation point ( x , y ) = ( 0 , 0 ), but unlike planar extensional flow, magnitudes of the velocity gradients in the x-direction (extensional) and y-direction (compressional) are unequal.

FIG. 2.

Flow field characterization for uniaxial and biaxial extensional flows (3D flows). (a) and (d) Visualization of 3D flow fields using simulations and experiments for uniaxial and biaxial extensional flow, respectively. Red lines indicate streamlines generated from finite-element simulations for uniaxial and biaxial extensional flow viewed in the x y plane using the identical flow geometries as in the experiments. Black lines denote particle trajectories determined from PTV experiments using 1  μm diameter particles in uniaxial and biaxial extensional flow, respectively, as viewed in the x y plane about z = 0. (b) and (e) Comparisons of velocity components v x and v y of particle trajectories as a function of position for x and y axes for simulations and experiments. The velocity components are determined in uniaxial and biaxial extensional flow with a 1.4 kPa pressure difference. (c) and (f) Comparison of d v x / d x and d v y / d y of PTV results of uniaxial and biaxial extensional flows as a function of applied pressures Δ p u and Δ p b, respectively. Dashed lines denote the idealized ratio for each flow field. Figure is shown in color in online version.

FIG. 2.

Flow field characterization for uniaxial and biaxial extensional flows (3D flows). (a) and (d) Visualization of 3D flow fields using simulations and experiments for uniaxial and biaxial extensional flow, respectively. Red lines indicate streamlines generated from finite-element simulations for uniaxial and biaxial extensional flow viewed in the x y plane using the identical flow geometries as in the experiments. Black lines denote particle trajectories determined from PTV experiments using 1  μm diameter particles in uniaxial and biaxial extensional flow, respectively, as viewed in the x y plane about z = 0. (b) and (e) Comparisons of velocity components v x and v y of particle trajectories as a function of position for x and y axes for simulations and experiments. The velocity components are determined in uniaxial and biaxial extensional flow with a 1.4 kPa pressure difference. (c) and (f) Comparison of d v x / d x and d v y / d y of PTV results of uniaxial and biaxial extensional flows as a function of applied pressures Δ p u and Δ p b, respectively. Dashed lines denote the idealized ratio for each flow field. Figure is shown in color in online version.

Close modal

Figure 2(b) shows particle velocities v x and v y as a function of distance from the stagnation point along the x- and y-axes. We determined the strain rates using a least squares regression to determine d v x / d x = ε u ˙. Here, the quantities d v x / d x and d v y / d y are determined from particle trajectories using PTV, and their ratio was determined to be 2.0 for the range of pressures examined here ( 1.4 Δ p u 6.9 kPa), which is in excellent agreement with the expected ratio of 2 for uniaxial extension as given by Eq. (1) [Fig. 2(c)].

Biaxial extensional flow is generated by applying equal positive pressures to the channels associated with the compressional z-axis ( p z and p + z) relative to the pressures applied to the four outlet channels ( p x, p + x, p y, p + y) [Figs. 2(d)2(f)]. We again characterized the flow field using the single-mode imaging setup [Fig. 1(c)]. Biaxial extensional flow is described by
v x = ε ˙ b x ,
(3a)
v y = ε ˙ b y ,
(3b)
v z = 2 ε ˙ b ,
(3c)
where ε ˙ b is the strain rate in the biaxial extensional flow. Similar to uniaxial extensional flow, particles exhibit nonzero v z velocities in biaxial extension. The number density of tracked particles is small near the stagnation point ( x , y ) = ( 0 , 0 ) because the trajectories of incoming particles from the compressional z-axis diverge from this point. The resulting flow field has radially diverging streamlines from the stagnation point when viewed in the x y plane (supplementary material, Movie 5 [61]). We determined d v x / d x and d v y / d y of the tracked particle trajectories and found that the ratios matched the expected value of unity within error [Fig. 2(f)], which validates the generation of biaxial extensional flow. Manipulation of the outlet pressures p x, p + x, p y, p + y can be used to tune this ratio to achieve arbitrary ellipsoidal extensional flows [26].

We next demonstrate trapping of single particles in 3D flow fields. Single colloidal particles (3.4  μm diameter) were trapped and confined near the stagnation point of uniaxial or biaxial extensional flow in low viscosity aqueous solution ( η = 1.0 mPa s). We began with particle trapping experiments in uniaxial extensional flow using the single-mode imaging setup (Fig. 3). A movie showing particle trapping in uniaxial extensional flow is provided in supplementary material (Movie 6 [61]). For these experiments, the detected particle center-of-mass position in the x-direction is compared to a desired setpoint position x 0, and the pressure difference ( Δ p x = p + x p x) is manipulated by the PID controller to adjust the stagnation point position to effectively confine particles in flow, similar to prior work on 2D flow trapping [Fig. 3(a)] [10,11]. Figure 3(c) shows the results from a long-time ( 10 min) trapping experiment, where the particle position along the extensional axis of uniaxial extensional flow is plotted as a function of time. The probability distribution of particle position along the x-direction during the 10-min trapping experiment is also plotted in Fig. 3(d), showing that the particle is effectively confined along the extensional axis. Trapped particles generally exhibit submicrometer standard deviation displacements from the mean along the unstable x-axis ( σ = 0.49 μm), which is significantly smaller than the bead diameter (3.4  μm). In addition, particles are passively trapped along the two orthogonal directions corresponding to the two compressional flow axes [Figs. 3(b), 3(e), and 3(f)]. For these experiments, a strain rate of ε ˙ u = 2.27 s 1 was used.

FIG. 3.

Particle trapping in uniaxial extensional flow. (a) Schematic showing particle trapping via manipulation of a stagnation point (denoted by cross) toward a setpoint position (open circle) in the x-direction (extensional axis). The offset Δ x between the particle position (filled circle) and setpoint position (open circle) is used to manipulate Δ p x to adjust the stagnation point position in the x-direction. (b) Schematic showing a particle passively trapped in the y-direction (compressional axis). (c) Time trace of x-position of the trapped particle. The inset schematic shows a trapped particle in uniaxial extensional flow, and the red arrows depict the x-axis. (d) Probability distribution of x-position from (c). The red line shows a Gaussian fit. (e) Time trace of y-position of the trapped particle. The inset schematic shows a trapped particle in uniaxial extensional flow, and the red arrows depict the y-axis. (f) Probability distribution of y-position from data shown in (e). The red line shows a Gaussian fit. (g) and (h) Power spectral densities of x- and y-position from data shown in (c) and (e) for the trapped particle. Red lines denote Lorentzian fits to Eq. (4). For these experiments, the uniaxial strain rate ε ˙ u = 2.27 s 1 and particles are trapped in low viscosity aqueous solutions ( η = 1.0 mPa s). Figure is shown in color in online version.

FIG. 3.

Particle trapping in uniaxial extensional flow. (a) Schematic showing particle trapping via manipulation of a stagnation point (denoted by cross) toward a setpoint position (open circle) in the x-direction (extensional axis). The offset Δ x between the particle position (filled circle) and setpoint position (open circle) is used to manipulate Δ p x to adjust the stagnation point position in the x-direction. (b) Schematic showing a particle passively trapped in the y-direction (compressional axis). (c) Time trace of x-position of the trapped particle. The inset schematic shows a trapped particle in uniaxial extensional flow, and the red arrows depict the x-axis. (d) Probability distribution of x-position from (c). The red line shows a Gaussian fit. (e) Time trace of y-position of the trapped particle. The inset schematic shows a trapped particle in uniaxial extensional flow, and the red arrows depict the y-axis. (f) Probability distribution of y-position from data shown in (e). The red line shows a Gaussian fit. (g) and (h) Power spectral densities of x- and y-position from data shown in (c) and (e) for the trapped particle. Red lines denote Lorentzian fits to Eq. (4). For these experiments, the uniaxial strain rate ε ˙ u = 2.27 s 1 and particles are trapped in low viscosity aqueous solutions ( η = 1.0 mPa s). Figure is shown in color in online version.

Close modal
Trap performance was characterized by determining the power spectral density (PSD) of the particle position fluctuations for a trapped particle [Figs. 3(g) and 3(h)] [9,46]. Briefly, the time traces of trapped particle positions were subdivided into N 8 bins with 50% overlap such that the length of each bin was an integer power of 2. A Hamming window was then applied to each bin, and the PSD was separately computed for each windowed bin, from which the average was determined, as shown in Figs. 3(g) and 3(h). The PSDs are well-described by a Lorentzian function
P ( f ) = k B T 2 π 2 ζ ( f c 2 + f 2 ) ,
(4)
where k B is Boltzmann’s constant, T is absolute temperature, ζ is the Stokes drag coefficient, f is the frequency, and f c is the corner frequency. Equation (4) is used to quantify trap performance by determining the trap stiffness, which is defined as κ = 2 π ζ f c. We note that Eq. (4) is derived on the basis where particles are passively trapped with a linear restoring force from a target setpoint, as is the case with optical traps. In the limit of small perturbations, the restoring force for a PID controller is approximated as a linear restoring force [11]. Using this approach, trap stiffness is determined as 1.4 × 10 4 pN/nm along the extensional axis [Fig. 3(g)], which is comparable to early generation hydrodynamic traps ( 1.2 × 10 4 pN/nm) [47]. In general, trap stiffness can be increased using an MPC scheme previously used in the Stokes trap for multiplexed particle manipulations [9]. Employing larger proportional controller gains generally resulted in particle oscillations about the setpoint, which manifests by the appearance of a peak in the PSD of particle position fluctuations [9,47].

Trapping particles in uniaxial extensional flow offers significant advantages over particle trapping in 2D planar extensional flow. In uniaxial extension, the second compressional flow axis orthogonal to the plane of imaging (i.e., along the z-direction) effectively pushes trapped particles into the plane of imaging, which prevents particles from diffusing out of the image plane in the z-direction. Planar extensional flow lacks this additional compressional flow axis orthogonal to the image plane, which necessitates active focusing during a trapping experiment due to particle diffusion in the z-direction. For trapping particles using 2D flows, particle diffusion inevitably results in particle motion in the z-direction and eventual approach to a confining boundary, which could affect the local strain rate in the 2D flow plane. Moreover, viscoelastic solutions with appreciable normal stress differences between the axis of extension and the neutral axis of planar extensional flow (e.g., τ x x τ z z) results in particle migration away from the center of the channel toward regions of lower strain rate [39,48,49]. Therefore, 3D flow fields greatly enhance long-time trap performance and stability compared to 2D flow traps due to passive particle confinement along the orthogonal z-direction, which is a major advantage of the 3D hydrodynamic trap that will be useful in achieving long-time observations of freely suspended active particles or objects in viscoelastic solutions [48,50].

We next trapped single colloidal particles in biaxial extensional flow using the single-mode imaging setup. A movie showing particle trapping in biaxial extensional flow is provided in supplementary material (Movie 7 [61]). Here, we utilized a similar feedback control scheme for uniaxial extension, but biaxial extensional flow has two extensional axes, which necessitates the use of two separate feedback controllers for the unstable x- and y-axes [Figs. 4(a) and 4(b)]. The use of two simultaneous controllers to manipulate Δ p x = p + x p x and Δ p y = p + y p y according to Δ x and Δ y offsets relies on the orthogonality or independence of setpoint changes between these two axes. Indeed, we observed that changes in Δ p x did not affect the stagnation point position in the y-direction and vice versa near the center of the cross-slot.

FIG. 4.

Particle trapping in biaxial extensional flow. (a) and (b) Schematics of particle trapping in biaxial extensional flow for the x- and y-axes, respectively. The stagnation point in biaxial extensional flow (denoted with a cross) is manipulated by modulating Δ p x and Δ p y based on the offsets between the particle position (filled circle) and the setpoint position (open circle), Δ x and Δ y, respectively. (c) Time trace of the x-position of the trapped bead. Inset schematic shows a trapped particle in biaxial extensional flow. The red arrows depict the x axis. (d) Probability distribution of x-position from data shown in (c). The red line is a Gaussian fit. (e) Time trace of y-position of the trapped bead. The inset schematic shows a trapped particle in the same biaxial extensional flow for the data shown in (c). The red arrows depict the y axis. (f) Probability distribution of y-position from data shown in (e). The red line is a Gaussian fit. (g) and (h) Power spectral densities of x- and y-position from (c) and (e) for the trapped particle. Red lines correspond to Lorentzian fits to Eq. (4). For these experiments, the biaxial strain rate ε ˙ b = 1.38 s 1, and particles are trapped in low viscosity aqueous solutions (1.0 mPa s). Figure is shown in color in online version.

FIG. 4.

Particle trapping in biaxial extensional flow. (a) and (b) Schematics of particle trapping in biaxial extensional flow for the x- and y-axes, respectively. The stagnation point in biaxial extensional flow (denoted with a cross) is manipulated by modulating Δ p x and Δ p y based on the offsets between the particle position (filled circle) and the setpoint position (open circle), Δ x and Δ y, respectively. (c) Time trace of the x-position of the trapped bead. Inset schematic shows a trapped particle in biaxial extensional flow. The red arrows depict the x axis. (d) Probability distribution of x-position from data shown in (c). The red line is a Gaussian fit. (e) Time trace of y-position of the trapped bead. The inset schematic shows a trapped particle in the same biaxial extensional flow for the data shown in (c). The red arrows depict the y axis. (f) Probability distribution of y-position from data shown in (e). The red line is a Gaussian fit. (g) and (h) Power spectral densities of x- and y-position from (c) and (e) for the trapped particle. Red lines correspond to Lorentzian fits to Eq. (4). For these experiments, the biaxial strain rate ε ˙ b = 1.38 s 1, and particles are trapped in low viscosity aqueous solutions (1.0 mPa s). Figure is shown in color in online version.

Close modal

Figures 4(c) and 4(e) show long-time traces ( 10 min) of particle position along the two extensional axes of biaxial extensional flow. In these experiments, a 3.4  μm diameter bead was confined in biaxial extensional flow at a strain rate ε ˙ b = 1.4 s 1. The probabilities of particle positions along the x- and y-directions (extensional direction) during the 10-min trapping experiment are plotted in Figs. 4(d) and 4(f). Similar to the case of uniaxial extensional flow, the trapped particle exhibits submicrometer standard deviation displacements in both the x and y-directions ( σ x = 0.57 and σ y = 0.50 μm), indicating robust trapping in two unstable axes over a period of 10 min. We further determined the PSD of particle position fluctuations in the x and y-directions and found that they are well-described by a Lorentzian given by Eq. (4) [Figs. 4(g) and 4(h)]. The experimentally determined trap stiffnesses of 1.0 × 10 4 and 1.9 × 10 4 pN/nm for the x- and y-axes are consistent with the trap stiffness determined from uniaxial extensional flow. To our knowledge, this work presents the first demonstration of particle trapping in biaxial extensional flows.

We next demonstrate 2D and 3D center-of-mass manipulation of particles using active flow control. In these experiments, the setpoint position is manipulated by independently modulating the pressure regulator pairs along the x, y, or z direction. We first demonstrate 2D particle manipulation in 3D flow fields using the single-mode imaging setup. Figures 5(a) and 5(b) show 2D center-of-mass particle trajectories in the x , y planes of uniaxial extension and biaxial extensional flow, respectively. In these experiments, a particle is passively localized along the compressional flow axes (i.e., y and z for uniaxial and z for biaxial extensional flow), and active feedback control is used to confine particles along the extensional axes, as previously described. To translate a particle along a user-defined trajectory in the x y plane, a target setpoint position is specified, and the feedback controller successively localizes the particles along the x and y directions to follow the setpoint. In this way, particles are manipulated along user-defined trajectories by successively specifying new trap setpoints along a trajectory in small increments [11]. We also demonstrate particle manipulation in the absence of an externally applied flow field [Fig. 5(c)]. By modulating Δ p x or Δ p y using the feedback controller, particles can be effectively manipulated even under conditions of zero-net imposed flow. Movies of particle manipulation in 2D using 3D flows are provided in supplementary material (Movies 1–3) [61].

FIG. 5.

Particle manipulation in 2D using 3D flow fields. Particle center-of-mass position is shown for manipulation in 3D flow fields. Closed particle trajectories were traced in 60 s along the clockwise direction. (a) Trajectory of the manipulated particle using uniaxial extensional flow at ε ˙ u = 2.27 s 1. (b) Trajectory of manipulated particle using biaxial extensional flow at ε ˙ b = 1.38 s 1. (c) Trajectory of manipulated particle using zero-net imposed external flow. Scale bars denote 50  μm. Figure is shown in color in online version.

FIG. 5.

Particle manipulation in 2D using 3D flow fields. Particle center-of-mass position is shown for manipulation in 3D flow fields. Closed particle trajectories were traced in 60 s along the clockwise direction. (a) Trajectory of the manipulated particle using uniaxial extensional flow at ε ˙ u = 2.27 s 1. (b) Trajectory of manipulated particle using biaxial extensional flow at ε ˙ b = 1.38 s 1. (c) Trajectory of manipulated particle using zero-net imposed external flow. Scale bars denote 50  μm. Figure is shown in color in online version.

Close modal

We further show full 3D positional control of particles using active flow control (Fig. 6). Here, we performed two different methods for localizing and manipulating particles in 3D. The particle z-position can be determined using the single-imaging mode coupled with an autofocusing scheme or by direct observation of the z-position using the dual-imaging configuration (Materials and methods section). For particle manipulation using the single-imaging mode, an autofocusing scheme is used to localize particles in the z-direction (supplementary material, Fig. 4 [61]) [39]. The average intensity I of a detected particle is compared against z-positions, and the direction of v z is determined from a computed gradient d I / d z such that the intensity of the particle is maximized to maintain focus. In brief, particle manipulation in the z-direction is achieved using an open-loop controller along the z-axis by applying a steady Δ p z offset after calibrating particle z locations with applied Δ p z. In the case of uniaxial extensional flow, this approach results in passive particle trapping in the z-direction. Supplementary material, Fig. 7 [61] shows particle trajectories following two different user-defined paths [supplementary material, Eq. (1) [61]] in the 3D space achieved by utilizing the single-mode autofocusing scheme in uniaxial extensional flow. In both cases, particle localization is generally more precise in the image plane ( x , y plane) compared to the z-direction due to the open-loop autofocusing method employed for determining particle focus in z. Particle localization along the y-compressional axis is expected to be comparable to the particle localization in the z-compressional axis in uniaxial extensional flow [Fig. 3(f)], where the particle experiences submicrometer standard-deviation of position fluctuations under a steady pressure setpoint. Therefore, we expect that particles are tightly confined in the z-direction, and the apparent variations in z-directions for the trajectories arise due to limited accuracy in determining the particle position using the fluorescence intensity-based autofocusing algorithm. Supplementary material, Fig. 8 [61] quantifies the particle deviation from the setpoint positions. Although the error in the particle position from the setpoint has a standard deviation of less than 1  μm in x and y, our results show that the autofocusing method resulted in a z-direction standard deviation of the position of >8  μm. Movies of particle manipulation in 3D using the single-imaging mode are provided in supplementary material (Movies 8 and 9)[61].

FIG. 6.

Particle manipulation in 3D using biaxial extensional flow at ε ˙ b = 0.59 s 1, respectively. (a) 3D particle trajectory during 3D manipulation using the orthogonal dual-imaging setup. Here, the particle moves along the z-direction when x and y are fixed to follow a prescribed trajectory. (b) 3D particle trajectory during center-of-mass manipulation in 3D using the dual-imaging mode. Here, x, y, and z are simultaneously changing during the particle translation experiment. Color bars denote the trajectory timepoints (units of seconds). Red traces indicate prescribed trajectories whereas colored traces denote observed particle coordinates. Figure is shown in color in online version.

FIG. 6.

Particle manipulation in 3D using biaxial extensional flow at ε ˙ b = 0.59 s 1, respectively. (a) 3D particle trajectory during 3D manipulation using the orthogonal dual-imaging setup. Here, the particle moves along the z-direction when x and y are fixed to follow a prescribed trajectory. (b) 3D particle trajectory during center-of-mass manipulation in 3D using the dual-imaging mode. Here, x, y, and z are simultaneously changing during the particle translation experiment. Color bars denote the trajectory timepoints (units of seconds). Red traces indicate prescribed trajectories whereas colored traces denote observed particle coordinates. Figure is shown in color in online version.

Close modal

To improve the localization accuracy for trapped particles in the z-direction, we implemented a dual-mode orthogonal-imaging setup to directly read out the particle z-position (Materials and methods section). Here, two orthogonal planes ( x y and y z planes) are imaged simultaneously to localize the particle center-of-mass in 3D (Fig. 6). In this configuration, the device position is fixed, and the x y and y z image plane objectives are continually translated to maintain the particle in focus in both x y and y z planes as the particle translates in 3D [supplementary material, Figs. 5 and 1(d) [61]]. Figure 6 shows two different center-of-mass trajectories of a particle along a user-defined path [supplementary material, Eq. (2) [61]] in the 3D space using biaxial extensional flow and the orthogonal dual-imaging setup. In the dual-imaging mode, the particle position in the z-direction is directly determined from center-of-mass localization of images in the y z plane, which leads to significant improvement in the accuracy of centroid localization of z-positions compared to the fluorescence intensity-based autofocusing scheme. In particular, trapped particles observed using the dual-imaging mode exhibit 1  μm standard deviation displacements from the setpoint position in the z-direction [supplementary material, Figs. 9(c) and 9(f)] [61]. Compared to the single-imaging mode [supplementary material, Figs. 8(c) and 8(f)] [61], the z-direction localization is significantly improved using the dual-imaging mode. Movies of particle manipulation in 3D using the dual-imaging mode are provided in supplementary material (Movies 10 and 11) [61].

As an initial demonstration of the new capabilities enabled by automated 3D flow control, we trapped and deformed liquid droplets in both uniaxial and biaxial extensional flows (Fig. 7). Water-in-oil droplets stabilized with surfactant were prepared off-device and introduced into the 3D cross-slot through the + z channel (Materials and methods section). In these experiments, the same droplet was deformed under both uniaxial and biaxial extensional flow. Strain rates were determined using the flow calibration data shown in supplementary material, Figs. 1 and 2 [61] while accounting for changes in solvent viscosity ( η o i l = 30 mPa s). Figures 7(a) and 7(b) show snapshots of a liquid droplet trapped in uniaxial extension at low ( ε ˙ u = 7 s 1) and high ( ε ˙ u = 34 s 1) strain rates, respectively. The droplet visibly stretches along the x-axis and deforms into a prolate spheroidal shape upon increasing the strain rate. Figures 7(c) and 7(d) show snapshots of the same droplet trapped in biaxial extensional flow at low ( ε ˙ b = 7 s 1) and high ( ε ˙ b = 43 s 1) strain rates. In biaxial extension, the radius of the trapped droplet increases in the x y plane increases as the droplet flattens along the z-axis and deforms into an oblate spheroid. In both uniaxial extension and biaxial extension, the droplet visibly deforms at higher strain rates, albeit with qualitatively different behavior.

FIG. 7.

Droplet deformation in uniaxial and biaxial extensional flows. (a) and (b) Images of a water-in-oil liquid droplet trapped in uniaxial extensional flow at ε ˙ u = 3 s 1 and at ε ˙ u = 34 s 1, respectively. (c) and (d) Images of a water-in-oil liquid droplet trapped in biaxial extensional flow at ε ˙ b = 5 s 1 and at ε ˙ b = 43 s 1, respectively. (e) Deformation parameter D of the droplet with increasing strain rates. The gray line denotes the linear fit of D in uniaxial flow, and the black line indicates the linear fit of D in biaxial flow. Error bars correspond to standard deviations.

FIG. 7.

Droplet deformation in uniaxial and biaxial extensional flows. (a) and (b) Images of a water-in-oil liquid droplet trapped in uniaxial extensional flow at ε ˙ u = 3 s 1 and at ε ˙ u = 34 s 1, respectively. (c) and (d) Images of a water-in-oil liquid droplet trapped in biaxial extensional flow at ε ˙ b = 5 s 1 and at ε ˙ b = 43 s 1, respectively. (e) Deformation parameter D of the droplet with increasing strain rates. The gray line denotes the linear fit of D in uniaxial flow, and the black line indicates the linear fit of D in biaxial flow. Error bars correspond to standard deviations.

Close modal
The droplet deformation parameter D was used to quantify the droplet shape in both uniaxial and biaxial extensional flows as follows [51]:
D = L B L + B ,
(5)
where L and B are the lengths of the major and minor axes of the ellipsoid, respectively. Although both L and B were directly measured from brightfield images in uniaxial extensional flow, only L could be measured from images in biaxial extension. Thus, the radius of the spherical droplet under zero flow conditions was estimated (supplementary material, Fig. 10 [61]), and volume conservation was used to calculate B, assuming that the droplet deformed into oblate spheroidal shape ( V L 2 B). The deformation parameters for a liquid droplet with radius R 0 = 84.0 μm at different uniaxial and biaxial strain rates are shown in Fig. 7(e). In all cases, D is observed to increase linearly with the strain rate in both flow regimes.
The extent of drop deformation is related to the magnitude of viscous drag forces relative to the surface tension forces, which is quantified by the capillary number C a. For extensional flow, C a is defined as
C a = η o i l R 0 ε ˙ γ .
(6)
Prior work based on numerical simulations has shown that in the limit of small C a, droplets exhibit equivalent deformation D in uniaxial and biaxial extensional flows such that D is related to C a through the relation [52,53]
D = 3 ( 16 + 19 λ ) 2 ( 16 + 16 λ ) C a ,
(7)
where λ is the ratio of the viscosity of the droplet to that of the surrounding liquid ( λ = η d r o p l e t / η o i l). Because the same droplet was used for both uniaxial and biaxial extensional flow experiments, D is expected to change nearly identically with increasing strain rate in both flow types. Indeed, dependence of D on ε ˙ is nearly identical for both uniaxial and biaxial extensional flows, as shown in Fig. 7(e). Although we did not have a priori knowledge of the surface tension of the droplet γ under the surfactant concentrations used, Eqs. (6) and (7) can be rearranged to show that d D / d ε ˙ 1 / γ. Experimental determination of d D / d ε ˙ [Fig. 7(e)] for the same droplet in uniaxial and biaxial extensional flows resulted in similar values of γ ( γ = 1.08 ± 0.03 and 1.17 ± 0.02 mN/m for the uniaxial and biaxial extensional flows, respectively). We note that Eqs. (6) and (7) assume Newtonian droplets with constant and uniform surface tension. However, in our case, the presence of surfactant is expected to cause some nonuniformity of surfactant concentrations along the droplet shape, where surfactants migrate toward the extensional axis [54]. Thus, the reported surface tension is an effective surface tension of a Newtonian droplet with uniform surfactant concentration. Utilizing these values of γ, we were able to determine the dimensionless flow strength C a. For the droplet deformation experiments shown here, C a < 0.1, indicating that the droplet was subjected to small deformations. At higher C a, droplets are expected to exhibit unique deformation dynamics in the two different flow fields, with interesting droplet breakup behavior including the possibility for equatorial tip streaming at large C a in biaxial extensional flow [52,55]. We anticipate that the 3D hydrodynamic trap will enable such studies in future work.
We further extended the technique to study viscoelastic materials by trapping and deforming lipid vesicles (GUVs) near the stagnation point of both uniaxial and biaxial extensional flows (Fig. 8). Vesicles were prepared off-device and introduced into the 3D cross-slot through the + z channel (Materials and methods section). Vesicle contours were extrapolated and revolved along the x-axis in uniaxial extensional flow to estimate the revolved volume and consequentially the reduced volume of the vesicles. For reduced volume estimation of the vesicles in biaxial extensional flow, the estimated length of the minor axis in the z-direction was considered before revolving the contours around the x-axis. (See the supplementary material [61].) Reduced volume ν is a measure of vesicle deflation and refers to the ratio of a vesicle’s volume V to the volume of an equivalent sphere with the same surface area A such that [44]
ν = 3 V 4 π R 3 ,
(8)
where R = A / 4 π is the vesicle’s equivalent radius based on the total surface area. Using the method described above, the reduced volumes of the vesicles in uniaxial extensional flow and biaxial extensional flow (Fig. 8) were estimated to be ν = 0.98 ± 0.06 and ν = 0.97 ± 0.06, respectively.
FIG. 8.

Vesicle deformation in uniaxial and biaxial extensional flows. (a) and (b) Images of nearly spherical DOPC vesicles in uniaxial extensional flow at ε ˙ u = 1.13 s 1 and at ε ˙ u = 4.17 s 1, respectively. (c) and (d) Images of DOPC vesicles trapped in biaxial extensional flow at ε ˙ b = 1.34 s 1 and at ε ˙ b = 2.24 s 1, respectively. (e) and (f) Steady-state deformation parameter D of two different vesicles as a function of capillary number C a in uniaxial extensional flow and biaxial extensional flow, respectively. The black lines denote linear fits of D. Error bars correspond to standard deviations in determining D.

FIG. 8.

Vesicle deformation in uniaxial and biaxial extensional flows. (a) and (b) Images of nearly spherical DOPC vesicles in uniaxial extensional flow at ε ˙ u = 1.13 s 1 and at ε ˙ u = 4.17 s 1, respectively. (c) and (d) Images of DOPC vesicles trapped in biaxial extensional flow at ε ˙ b = 1.34 s 1 and at ε ˙ b = 2.24 s 1, respectively. (e) and (f) Steady-state deformation parameter D of two different vesicles as a function of capillary number C a in uniaxial extensional flow and biaxial extensional flow, respectively. The black lines denote linear fits of D. Error bars correspond to standard deviations in determining D.

Close modal
For vesicle dynamics in extensional flow, a stability boundary was predicted [56] and experimentally verified [44] based on the relationship between reduced volume ν, viscosity contrast λ, and capillary number C a, defined as [56]
C a = η ε ˙ R 3 κ b ,
(9)
where η is the viscosity of the fluid surrounding the vesicle and κ b is the membrane’s bending modulus. The viscosity contrast λ is the ratio of the fluid viscosities between the interior and the exterior regions of a vesicle ( λ = η i n / η o u t). In this study, we observed the nonequilibrium dynamics of vesicles in both uniaxial and biaxial extensional flows over a range of capillary number C a for a viscosity contrast λ = 1.

The scope of the current study was to demonstrate the 3D manipulation and deformation of lipid vesicles in 3D flows. Here, we focused on nondeflated, nearly spherical vesicles ( ν 1) that generally show significantly less deformation in 2D planar extensional flow compared to deflated vesicles ( ν 0.75) [44]. For the bending modulus, κ b was taken from previously determined values for DOPC lipid vesicles with 0.12 mol. % DOPE-Rh and 100 mM sucrose at T = 24 °C ( κ b = 9.17 × 10 20 J) [44]. Although the DOPC vesicles in this work are labeled with Cy5 dye instead of rhodamine dye, prior work has shown that the small amount of DOPE-Cy5 labeled lipid does not significantly alter the bending modulus of the membrane compared to pure DOPC vesicles [57].

Nearly spherical vesicles were trapped in uniaxial and biaxial extensional flows, and the deformation was quantified as a function of C a. Here, the deformation parameter D [Eq. (5)] was used to quantify vesicle deformation in both uniaxial and biaxial extensional flows. Figures 8(a) and 8(b) show snapshots of a vesicle trapped in uniaxial extension at ε ˙ u = 1.13 s 1 and ε ˙ u = 4.17 s 1, respectively. As the strain rate increases, the vesicle visibly elongates along the x-axis and undergoes deformation into a prolate spheroidal shape. Figures 8(c) and 8(d) show snapshots of another vesicle trapped in biaxial extension at ε ˙ b = 1.34 s 1 and ε ˙ b = 2.24 s 1, respectively. Movies of vesicles stretching in uniaxial and biaxial extensional flows are provided in supplementary material (Movies 14 and 15) [61]. The radius of each vesicle was determined from the images using volume conservation assuming that the vesicles deformed into a prolate spheroidal shape and an oblate spheroidal shape in the uniaxial and biaxial extensional flow, respectively ( V L 2 B). The steady-state deformation parameters for two nearly spherical vesicles with R o = 15.3 and R o = 24.6 μm as a function of C a are shown in Figs. 8(e) and 8(f). In all cases, D was observed to increase linearly with capillary number in both flow regimes.

In this work, we demonstrate trapping and 2D and 3D manipulation of small colloidal particles in 3D uniaxial and biaxial extensional flows. We characterized the 3D flow fields under the steady pressure-driven flow generated in a 6-arm cross-slot geometry. Quantitative analysis of the particle velocities in the x and y axes show that the flow fields match the expected velocity profiles for uniaxial and biaxial extensional flows. Active feedback control is used to confine and manipulate particles along the unstable (extensional) flow axes in uniaxial and biaxial extensional flows, and stable particle trapping is demonstrated for micrometer-sized particles for long trapping durations ( 10 min). In all cases, characterization of the trap stiffnesses shows that the trapping performance of the 3D flow device is comparable to prior 2D flow trapping devices. For the first time, we demonstrate successful flow-based manipulation of the particle center-of-mass position in 2D and 3D by moving particles along user-defined trajectories coupled with active feedback control. We further demonstrate an application of particle trapping in 3D flow fields by deforming liquid droplets and GUVs in both uniaxial and biaxial extensional flows and find nearly identical linear deformation behavior for a given microdroplet under both flow fields.

Overall, the results presented in this work extend the capabilities of flow-based trapping to manipulate particles in three dimensions. The 3D hydrodynamic trap presented here can be further improved with more advanced numerical models and control mechanisms to study the dynamics of soft matter in 3D flows [58]. The studies of vesicle and droplet deformation show initial examples of how 3D cross-slot flow devices will enable new studies of soft materials in nonequilibrium flow, including deformation dynamics that cannot be observed in conventional 2D microfluidic devices [55]. Moreover, our results demonstrate that using this technique, we can successfully trap and manipulate GUVs in both uniaxial and biaxial extensional flows. This represents important implications for the future study of biological membranes and other viscoelastic materials in nonequilibrium flow. In addition, the 3D printing technique used here enables facile construction of 3D geometries consisting of an arbitrary number of channels, which will enable control over more state variables in addition to the center-of-mass position of one particle as shown in this work (e.g., multiplexed particle manipulation in 3D) [9]. To our knowledge, this is the first reported method for trapping single particles in biaxial extensional flow, which is rarely studied in microfluidic systems. We anticipate that the 3D hydrodynamic trap will open new avenues in studying viscoelastic materials that may respond to more than one extensional axis, such as ring polymers which are expected to open along the extensional plane of biaxial extension due to their unique closed topology [59,60]. Future work will enable detailed investigations of the dynamics of soft materials in full 3D flow fields that more closely resemble common industrial processing conditions.

This work was supported by the National Science Foundation under Grant No. CBET-2030537, the National Institutes of Health under Grant No. 1-R01-GM143723-01A1 (E.F. and C.M.S.), the Joint Center for Energy Storage Research (JCESR), an Energy Innovation Hub funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (H.V.N. and C.M.S.), a Beckman Institute Postdoctoral Fellowship (M.I.J.), a Mistletoe Fellowship (M.I.J.), and a Sandia Academic Alliance fellowship (M.Q.T. and C.M.S.).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
Ashkin
,
A.
,
J. M.
Dziedzic
,
J. E.
Bjorkholm
, and
S.
Chu
, “
Observation of a single-beam gradient force optical trap for dielectric particles
,”
Opt. Lett.
11
,
288
290
(
1986
).
2.
Strick
,
T. R.
,
J.-F.
Allemand
,
D.
Bensimon
,
A.
Bensimon
, and
V.
Croquette
, “
The elasticity of a single supercoiled DNA molecule
,”
Science
271
,
1835
1837
(
1996
).
3.
Hertz
,
H.
, “
Standing-wave acoustic trap for nonintrusive positioning of microparticles
,”
J. Appl. Phys.
78
,
4845
4849
(
1995
).
4.
Cohen
,
A. E.
, and
W.
Moerner
, “
Method for trapping and manipulating nanoscale objects in solution
,”
Appl. Phys. Lett.
86
,
093109
(
2005
).
5.
Bentley
,
B.
, and
L.
Leal
, “
A computer-controlled four-roll mill for investigations of particle and drop dynamics in two-dimensional linear shear flows
,”
J. Fluid Mech.
167
,
219
240
(
1986
).
6.
Schroeder
,
C. M.
, “
Single polymer dynamics for molecular rheology
,”
J. Rheol.
62
,
371
403
(
2018
).
7.
Kumar
,
D.
,
A.
Shenoy
,
S.
Li
, and
C. M.
Schroeder
, “
Orientation control and nonlinear trajectory tracking of colloidal particles using microfluidics
,”
Phys. Rev. Fluids
4
,
114203
(
2019
).
8.
Kumar
,
D.
,
A.
Shenoy
,
J.
Deutsch
, and
C. M.
Schroeder
, “
Automation and flow control for particle manipulation
,”
Curr. Opin. Chem. Eng.
29
,
1
8
(
2020
).
9.
Shenoy
,
A.
,
C. V.
Rao
, and
C. M.
Schroeder
, “
Stokes trap for multiplexed particle manipulation and assembly using fluidics
,”
Proc. Natl. Acad. Sci. U.S.A.
113
,
3976
3981
(
2016
).
10.
Tanyeri
,
M.
,
M.
Ranka
,
N.
Sittipolkul
, and
C. M.
Schroeder
, “
A microfluidic-based hydrodynamic trap: Design and implementation
,”
Lab Chip
11
,
1786
1794
(
2011
).
11.
Tanyeri
,
M.
, and
C. M.
Schroeder
, “
Manipulation and confinement of single particles using fluid flow
,”
Nano Lett.
13
,
2357
2364
(
2013
).
12.
Brimmo
,
A. T.
, and
M. A.
Qasaimeh
, “
Stagnation point flows in analytical chemistry and life sciences
,”
RSC Adv.
7
,
51206
51232
(
2017
).
13.
Shenoy
,
A.
,
M.
Tanyeri
, and
C. M.
Schroeder
, “
Characterizing the performance of the hydrodynamic trap using a control-based approach
,”
Microfluid. Nanofluid.
8
,
1055
1066
(
2015
).
14.
Love
,
J. C.
,
J. R.
Anderson
, and
G. M.
Whitesides
, “
Fabrication of three-dimensional microfluidic systems by soft lithography
,”
MRS Bull.
26
,
523
528
(
2001
).
15.
Anna
,
S. L.
,
G. H.
McKinley
,
D. A.
Nguyen
,
T.
Sridhar
,
S. J.
Muller
,
J.
Huang
, and
D. F.
James
, “
An interlaboratory comparison of measurements from filament-stretching rheometers using common test fluids
,”
J. Rheol.
45
,
83
114
(
2001
).
16.
Bach
,
A.
,
H. K.
Rasmussen
, and
O.
Hassager
, “
Extensional viscosity for polymer melts measured in the filament stretching rheometer
,”
J. Rheol.
47
,
429
441
(
2003
).
17.
Chatraei
,
S.
,
C. W.
Macosko
, and
H.
Winter
, “
Lubricated squeezing flow: A new biaxial extensional rheometer
,”
J. Rheol.
25
,
433
443
(
1981
).
18.
Engmann
,
J.
,
C.
Servais
, and
A. S.
Burbidge
, “
Squeeze flow theory and applications to rheometry: A review
,”
J. Non-Newtonian Fluid Mech.
132
,
1
27
(
2005
).
19.
Teo
,
W. E.
, and
S.
Ramakrishna
, “
A review on electrospinning design and nanofibre assemblies
,”
Nanotechnology
17
,
R89
R106
(
2006
).
20.
Persano
,
L.
,
A.
Camposeo
,
C.
Tekmen
, and
D.
Pisignano
, “
Industrial upscaling of electrospinning and applications of polymer nanofibers: A review
,”
Macromol. Mater. Eng.
298
,
504
520
(
2013
).
21.
Barborik
,
T.
, and
M.
Zatloukal
, “
Steady-state modeling of extrusion cast film process, neck-in phenomenon, and related experimental research: A review
,”
Phys. Fluids
32
,
061302
(
2020
).
22.
Choi
,
K.-J.
,
J. E.
Spruiell
, and
J. L.
White
, “
Orientation and morphology of high-density polyethylene film produced by the tubular blowing method and its relationship to process conditions
,”
J. Polym. Sci.: Polym. Phys. Ed.
20
,
27
47
(
1982
).
23.
van Berkel
,
J. G.
,
N.
Guigo
,
J. J.
Kolstad
, and
N.
Sbirrazzuoli
, “
Biaxial orientation of poly (ethylene 2, 5-furandicarboxylate): An explorative study
,”
Macromol. Mater. Eng.
303
,
1700507
(
2018
).
24.
Nishioka
,
A.
,
T.
Takahashi
,
Y.
Masubuchi
,
J.-I.
Takimoto
, and
K.
Koyama
, “
Description of uniaxial, biaxial, and planar elongational viscosities of polystyrene melt by the K-BKZ model
,”
J. Non-Newtonian Fluid Mech.
89
,
287
301
(
2000
).
25.
Wagner
,
M.
, and
J.
Schaeffer
, “
Nonlinear strain measures for general biaxial extension of polymer melts
,”
J. Rheol.
36
,
1
26
(
1992
).
26.
Haward
,
S.
,
C.
Hopkins
,
K.
Toda-Peters
, and
A. Q.
Shen
, “
Microfluidic analog of an opposed-jets device
,”
Appl. Phys. Lett.
114
,
223701
(
2019
).
27.
Fuller
,
G. G.
,
C. A.
Cathey
,
B.
Hubbard
, and
B. E.
Zebrowski
, “
Extensional viscosity measurements for low-viscosity fluids
,”
J. Rheol.
3
,
235
249
(
1987
).
28.
Cathey
,
C. A.
, and
G. G.
Fuller
, “
Uniaxial and biaxial extensional viscosity measurements of dilute and semi-dilute solutions of rigid rod polymers
,”
J. Non-Newtonian Fluid Mech.
30
,
303
316
(
1988
).
29.
Hermansky
,
C. A.
, and
D. V.
Boger
, “
Opposing-jet viscometry of fluids with viscosity approaching that of water
,”
J. Non-Newtonian Fluid Mech.
56
,
1
14
(
1995
).
30.
Gonzalez
,
J.
, and
B.
Liu
, “
Symmetry-based nonperturbative micromanipulation in a three-dimensional microfluidic device
,”
Phys. Rev. Fluids
5
,
044202
(
2020
).
31.
Afonso
,
A.
,
M.
Alves
, and
F.
Pinho
, “
Purely elastic instabilities in three-dimensional cross-slot geometries
,”
J. Non-Newtonian Fluid Mech.
165
,
743
751
(
2010
).
32.
Cruz
,
F.
,
R.
Poole
,
A.
Afonso
,
F.
Pinho
,
P.
Oliveira
, and
M.
Alves
, “
Influence of channel aspect ratio on the onset of purely-elastic flow instabilities in three-dimensional planar cross-slots
,”
J. Non-Newtonian Fluid Mech.
227
,
65
79
(
2016
).
33.
Zografos
,
K.
,
N.
Burshtein
,
A. Q.
Shen
,
S. J.
Haward
, and
R. J.
Poole
, “
Elastic modifications of an inertial instability in a 3D cross-slot
,”
J. Non-Newtonian Fluid Mech.
262
,
12
24
(
2018
).
34.
Burshtein
,
N.
,
S. T.
Chan
,
K.
Toda-Peters
,
A. Q.
Shen
, and
S. J.
Haward
, “
3D-printed glass microfluidics for fluid dynamics and rheology
,”
Curr. Opin. Colloid Interface Sci.
43
,
1
14
(
2019
).
35.
Sun
,
H.
,
Y.
Jia
,
Y.
Dong
,
D.
Dong
, and
Z.
J.
, “
Combining additive manufacturing with microfluidics: An emerging method for developing novel organs-on-chips
,”
Curr. Opin. Chem. Eng.
28
,
1
9
(
2020
).
36.
Lee
,
W.
,
D.
Kwon
,
W.
Choi
,
G. Y.
Jung
,
A. K.
Au
,
A.
Folch
, and
S.
Jeon
, “
3D-printed microfluidic device for the detection of pathogenic bacteria using size-based separation in helical channel with trapezoid cross-section
,”
Sci. Rep.
5
,
1
7
(
2015
).
37.
Spivey
,
E. C.
,
B.
Xhemalce
,
J. B.
Shear
, and
I. J.
Finkelstein
, “
3D-printed microfluidic microdissector for high-throughput studies of cellular aging
,”
Anal. Chem.
86
,
7406
7412
(
2014
).
38.
Anthony
,
S.
,
L.
Zhang
, and
S.
Granick
, “
Methods to track single-molecule trajectories
,”
Langmuir
22
,
5266
5272
(
2006
).
39.
Hsiao
,
K.-W.
,
J.
Dinic
,
Y.
Ren
,
V.
Sharma
, and
C. M.
Schroeder
, “
Passive non-linear microrheology for determining extensional viscosity
,”
Phys. Fluids
29
,
121603
(
2017
).
40.
Shum
,
H.
,
A.
Sauret
,
A.
Fernandez-Nieves
,
H.
Stone
, and
D.
Weitz
, “
Corrugated interfaces in multiphase core-annular flow
,”
Phys. Fluids
22
,
082002
(
2010
).
41.
Tsai
,
S.
,
J.
Wexler
,
J.
Wan
, and
H.
Stone
, “
Microfluidic ultralow interfacial tensiometry with magnetic particles
,”
Lab Chip
13
,
119
125
(
2013
).
42.
Motagamwala
,
A. H.
, A microfluidic, extensional flow device for manipulating soft particles, Ph.D. thesis, Department of Chemical Engineering and Applied Chemistry, University of Toronto, 2013.
43.
Angelova
,
M.
,
S.
Soléau
,
P.
Méléard
,
F.
Faucon
, and
P.
Bothorel
, Preparation of giant vesicles by external AC electric fields. Kinetics and applications, in
Trends Colloid Interface Science VI
(Steinkopff, Darmstadt, 1992), pp. 127–131.
44.
Kumar
,
D.
,
C. M.
Richter
, and
C. M.
Schroeder
, “
Conformational dynamics and phase behavior of lipid vesicles in a precisely controlled extensional flow
,”
Soft Matter
16
,
337
347
(
2020
).
45.
Kumar
,
D.
, and
C. M.
Schroeder
, “
Nonlinear transient and steady-state stretching of deflated vesicles in flow
,”
Langmuir
37
,
13976
13984
(
2021
).
46.
Lansdorp
,
B. M.
, and
O. A.
Saleh
, “
Power spectrum and allan variance methods for calibrating single-molecule video-tracking instruments
,”
Rev. Sci. Instrum.
83
,
025115
(
2012
).
47.
Tanyeri
,
M.
,
E. M.
Johnson-Chavarria
, and
C. M.
Schroeder
, “
Hydrodynamic trap for single particles and cells
,”
Appl. Phys. Lett.
96
,
224101
(
2010
).
48.
Hsiao
,
K.-W.
,
C.
Sasmal
,
J.
Ravi Prakash
, and
C. M.
Schroeder
, “
Direct observation of dna dynamics in semidilute solutions in extensional flow
,”
J. Rheol.
61
,
151
167
(
2017
).
49.
Patel
,
S. F.
,
C. D.
Young
,
C. E.
Sing
, and
C. M.
Schroeder
, “
Nonmonotonic dependence of comb polymer relaxation on branch density in semidilute solutions of linear polymers
,”
Phys. Rev. Fluids
5
,
121301
(
2020
).
50.
Johnson-Chavarria
,
E. M.
,
U.
Agrawal
,
M.
Tanyeri
,
T. E.
Kuhlman
, and
C. M.
Schroeder
, “
Automated single cell microbioreactor for monitoring intracellular dynamics and cell growth in free solution
,”
Lab Chip
14
,
2688
2697
(
2014
).
51.
Taylor
,
G. I.
, “
The formation of emulsions in definable fields of flow
,”
Proc. R. Soc. London, Ser. A
146
,
501
523
(
1934
).
52.
Stone
,
H.
, and
L.
Leal
, “
A note concerning drop deformation and breakup in biaxial extensional flows at low Reynolds numbers
,”
J. Colloid Interface Sci.
133
,
340
347
(
1989
).
53.
Favelukis
,
M.
, “
A drop in uniaxial and biaxial nonlinear extensional flows
,”
Phys. Fluids
29
,
087102
(
2017
).
54.
Suryo
,
R.
, and
O. A.
Basaran
, “
Tip streaming from a liquid drop forming from a tube in a co-flowing outer fluid
,”
Phys. Fluids
18
,
082102
(
2006
).
55.
Brosseau
,
Q.
, and
P.
Vlahovska
, “
Streaming from the equator of a drop in an external electric field
,”
Phys. Rev. Lett.
119
,
034501
(
2017
).
56.
Narsimhan
,
V.
,
A. P.
Spann
, and
E. S.
Shaqfeh
, “
The mechanism of shape instability for a vesicle in extensional flow
,”
J. Fluid Mech.
750
,
144
190
(
2014
).
57.
Bouvrais
,
H.
,
T.
Pott
,
L. A.
Bagatolli
,
J. H.
Ipsen
, and
P.
Méléard
, “
Impact of membrane-anchored fluorescent probes on the mechanical properties of lipid bilayers
,”
Biochim. Biophys. Acta
1798
,
1333
1337
(
2010
).
58.
Fang
,
W.-Z.
,
T.
Xiong
,
O. S.
Pak
, and
L.
Zhu
, “
Data-driven intelligent manipulation of particles in microfluidics
,”
Adv. Sci.
10
,
2205382
(
2023
).
59.
Tu
,
M.
,
M.
Lee
,
R. M.
Robertson-Anderson
, and
C. M.
Schroeder
, “
Direct observation of ring polymer dynamics in the flow-gradient plane of shear flow
,”
Macromolecules
53
,
9406
9419
(
2020
).
60.
Murashima
,
T.
,
K.
Hagita
, and
T.
Kawakatsu
, “
Viscosity overshoot in biaxial elongational flow: Coarse-grained molecular dynamics simulation of ring–linear polymer mixtures
,”
Macromolecules
54
,
7210
7225
(
2021
).
61.
See supplementary material online for details on flow field characterization, z-position localization and 3D trajectories, dual-imaging setup, feedback control algorithm, live image processing, and droplet and vesicle deformation experiments.

Supplementary Material