Viscous flow through microfabricated axisymmetric contraction/expansion geometries

We employ a state-of-the-art microfabrication technique (selective laser-induced etching) to fabricate a set of axisymmetric microfluidic geometries featuring a 4:1 contraction followed by a 1:4 downstream expansion in the radial dimension. Three devices are fabricated: the first has a sudden contraction followed by a sudden expansion, the second features hyperbolic contraction and expansion profiles, and the third has a numerically optimized contraction/expansion profile intended to provide a constant extensional/compressional rate along the axis. We use micro-particle image velocimetry to study the creeping flow of a Newtonian fluid through the three devices and we compare the obtained velocity profiles with finite-volume numerical predictions, with good agreement. This work demonstrates the capability of this new microfabrication technique for producing accurate non-planar microfluidic geometries with complex shapes and with sufficient clarity for optical probes. The axisymmetric microfluidic geometries examined have potential to be used for the study of the extensional properties and non-linear dynamics of viscoelastic flows, and to investigate the transport and deformation dynamics of bubbles, drops, cells, and fibers.


Introduction
Microfluidics involves fluid flows at small length scales, typically ∼ 100 μ m (Stone and Kim 2001). Over the last two decades, the fabrication of microfluidic devices has been dominated by soft lithography techniques, which enable precise fabrication of channels with a planar rectangular cross-section, with complex variations of the geometry only possible in two dimensions (2D) (McDonald and Whitesides 2002). However, in many cases, a 2D approximation of a problem of interest is not ideal and researchers would prefer, if possible, to retain the use of the third dimension in the design of their model microfluidic systems. A prime example is for the in vitro study of hemodynamics, which requires model blood vessels with a circular cross-section (Doutel et al. 2015;Kang et al. 2018). In addition, many processes of interest, for example spraying, inkjet printing, and fiber spinning, involve fluid flows through circular ducts and nozzles or spinnerets. In all these processes, fluid elements are subjected to extensional or compressional velocity gradients at locations where the cross-sectional area of a duct varies. A feature of the small length scales of microfluidics is the ability to reach high deformation rates ( ∼ −1 ) for relatively weak (even negligible) inertial forces (since Reynolds numbers Re ∼ ) (Squires and Quake 2005). At the microscale, fluids with even weakly non-Newtonian properties such as viscoelasticity (exhibited by, e.g., blood, paints, and printing inks) can display strongly non-linear behavior due to any extensional kinematics in the flow field. This occurs due to the non-linear increase in elastic tensile stresses (or the extensional viscosity of the fluid) in regions of high strain rates and high fluid strains, resulting in feedback on the flow field. In fact, the extensional viscosity of a viscoelastic fluid often dominates the macroscopic flow response at small length scales, and is, thus, a material function of fundamental importance in a wide range of processes (Galindo-Rosales et al. 2013;Keshavarz and McKinley 2016;Haward 2016). The extensional viscosity of even simple Newtonian fluids is theoretically predicted to depend on the mode of extensional deformation: uniaxial, planar, or biaxial (Petrie 2006). This difference is likely to be greatly pronounced for viscoelastic fluids, so the development of axisymmetric microdevices for extensional rheometry of complex fluids is a worthwhile endeavor.
Among the various methods proposed for the measurement of extensional properties (or extensional rheology) of viscoelastic fluids, perhaps, the simplest is to use the extensional flow generated by a change in the cross-sectional area of the flow field. At the macroscale, there is a long history of using planar and axisymmetric contraction/expansion geometries for examining the extensional rheology of viscous and highly non-Newtonian fluids (e.g., Cogswell 1978;Evans and Walters 1986;Walters and Webster 1982;James et al. 1990;McKinley et al. 1991;Rothstein and McKinley 1999;Nigen and Walters 2002). These studies show that the response of a given fluid can be significantly different in an axisymmetric compared to a planar geometry. However, to date, microscale studies (more suited to low viscosity and weakly non-Newtonian fluids) have been restricted to simple planar geometries. Initially, planar abrupt microfluidic contractions with a contraction ratio (the ratio of upstream to downstream channel cross-sectional area) C r = 16 were employed by Rodd et al. (2005Rodd et al. ( , 2007, while Li et al. (2011a, b) extended the study to C r = 4 and 8. The significance of C r is that it defines the logarithmic (Hencky) strain on fluid elements traversing the contraction H = ln(C r ) . These early studies with polymeric fluids largely focussed on the vortex growth and non-linear dynamics upstream of the contraction plane, but also demonstrated the potential of the geometry for extensional rheometry purposes. However, a major drawback of an abrupt contraction is the non-uniform velocity gradient across the contraction. This is, in principle, resolved by making a hyperbolic converging channel wall profile, as described by Oliveira et al. (2007). These authors performed Newtonian flow experiments and numerical simulations to examine the homogeneity of the extensional rate along the flow axis in converging planar microgeometries with various hyperbolic profiles. A planar hyperbolic contraction/expansion geometry has since been commercialized as a high deformation rate extensional viscometer rheometer on-a-chip device (the e-VROC TM , Rheosense Inc., CA, see e.g., Pipe et al. 2008;Ober et al. 2013;Keshavarz and McKinley 2016;Garcia and Saraji 2019;Zografos et al. 2020). Similar microfluidic hyperbolic converging channels have also been employed to simulate blood flows in stenoses (Sousa et al. 2011), to measure the deformation of white blood cells (Rodrigues et al. 2015) and single DNA molecules (Larson et al. 2006), to study the deformation and breakup of droplets (Mulligan and Rothstein 2011), and to examine elastic turbulence of viscoelastic fluids (Groisman and Steinberg 2000;Ekanem et al. 2020). Some of the most recent developments in this area have involved "fine-tuning" the hyperbolic wall profile of the planar converging/ diverging geometry by an iterative numerical optimization procedure to further improve the homogeneity of the flow field (Zografos et al. 2016). These planar optimized devices have already found applications in the study of the buckling dynamics of individual suspended actin filaments and for the mechanical characterization of single DNA molecules and protein aggregates (Chakrabarti et al. 2020;Liu et al. 2020). We note that it is not only the extensional viscosity of a complex fluid that is expected to depend on the applied mode of elongation. Also, since the compression/elongation in uni/biaxial extension is along two orthogonal axes (with elongation/compression along the third direction), fundamental differences can be expected in the morphologies and dynamics observed for soft objects in axisymmetric compared with planar extension (where one direction is neutral).
In this work, we aim to extend the study of microfluidic contraction/expansion flows to the axisymmetric case, which more closely represent the kinematics found in real biological systems and flow loops in general. For this, we employ a novel microfabrication technique to etch channels of circular cross-section in fused-silica glass, with high precision. We begin with the relatively simple abrupt axisymmetric contraction/expansion, before examining a slowly converging/diverging hyperbolic-shaped channel. Finally, we fabricate and test a numerically optimized geometry designed to improve the homogeneity of the extensional flow field. We examine viscous Newtonian flow in the channels using complementary experimental and numerical methods. While this is obviously a first step, nevertheless it provides a valuable characterization of the flow field for assessing the utility of such devices for purposes such as the study of the non-linear dynamics and extensional rheometry of viscoelastic fluids or the dynamics of deformable objects. Furthermore, our efforts here can lead to further advances in the fabrication of three-dimensional (3D) microfluidic devices and to more complex 3D numerical optimizations of, e.g., junctions and flow-focusing devices for applications such as droplet and fiber formation (Pimenta et al. 2018;Zografos et al. 2019). One of our ultimate goals in this line of research is the optimization of the 3D 6-arm cross-slot device, capable of generating uniaxial and general biaxial extensional flows with a free stagnation point and unlimited Hencky strain (Afonso et al. 2010;Haward et al. 2012Haward et al. , 2019.

Axisymmetric microfluidic devices
General description In this work, we employ three microfabricated axisymmetric geometries featuring a contraction and subsequent expansion in the cross-sectional area. The three geometries (abrupt, hyperbolic, and optimized) are shown in schematic form in Fig. 1. We define the coordinate system (r, , z) , where r is the radial position, is the azimuthal angle, z is the axial (streamwise) position, and the origin is located at the center of the contraction throat. For all three channels, the characteristic radius in the up/downstream in/outlets is R in and the radius decreases to a minimum value of R c = 0.25R in . The contraction ratio is thus C r = (R in ∕R c ) 2 = 16 and the Hencky strain imposed in all three channels is H = ln 16 ≈ 2.77.
The surface wall profile of the abrupt contraction/expansion is described by: that is, the narrow constriction with constant radius R c spans a length L c = 9R in .
The hyperbolic contraction/expansion geometry converges/ diverges over the same distance of L c = 9R in and its surface wall profile can be described as: Since the channel cross-sectional area A(z) = R(z) 2 , and for Poiseuille flow in a circular duct, the streamwise flow velocity at r = 0 is v z (z) = 2Q∕A(z) (where Q = R 2 in U in is the imposed volumetric flow rate and U in is the average flow velocity across the inlet and outlet channels), it can be shown that the extensional rate along the flow axis of the channel ̇(z) = v z (z)∕ z is (nominally) given by: Fig. 1 Schematic diagrams of the three axisymmetric contraction/expansion geometries, showing the coordinate system and labels for the principal dimensions Note that the nominal extension rate profile given by Eq. 3 is not expected to be entirely accurate, since it assumes that the flow fully develops instantaneously at each axial location and that there is no formation of Moffatt-like vortices (Moffatt 1964) in the salient corners upstream of the contraction and downstream of the expansion (Oliveira et al. 2007).
Our approach to produce a closer approximation to the nominal extension rate profile described by Eq. 3 is by the modification of the simple hyperbolic shape via numerical optimization of the channel design. The optimized channel shape is determined by the numerical procedure described in Sect. 3.1 performed over a length L opt = 13R in (i.e., spanning the range −6.5R in ≤ z ≤ 6.5R in ). The resulting wall profile is complex and cannot be described by a simple piecewise function, although (as can be seen from Fig. 1) the profile collapses towards the hyperbolic shape close to z = 0 . The modification of the optimized shape from the hyperbolic is mostly characterized by the bulges that appear for |z| > 4.5R in . Comparable geometric modifications were found to result from optimizations performed with planar contraction/expansion geometries (Zografos et al. 2016).
Fabrication The microfluidic devices are fabricated in fused-silica glass by the subtractive 3D-printing method of selective laser-induced etching (SLE) (Gottmann et al. 2012;Meineke et al. 2016;Burshtein et al. 2019). SLE involves the modification of a volume of a fused-silica substrate, according to a specified 3D microchannel design, by exposure to a femtosecond-pulsed laser. Subsequently, the modified volume of material is selectively removed (approximately 1000 times faster than the unmodified material) by a chemical etchant (potassium hydroxide at 80 • C). The resulting microchannels are contained within a monolithic piece of highly transparent, rigid, and chemically resistant material, capable of withstanding high pressures without deforming significantly under flow. For the channel fabrication, we employed the commercial LightFab SLE printer (LightFab GmbH, Germany), which accommodates fused-silica substrates of the size of standard laboratory glass slides ( ≈ 25 × 75 mm), and with a maximum thickness of 5 mm along the direction of the light path. The instrument allows the fabrication of microchannels of cross-sectional dimension down to O(10 μ m) with arbitrary 3D shapes (no limit to the angle of undercuts), and a precision and root-mean-square surface roughness of O(1 μm). We note that slight channel tapering can occur for very long channels (e.g., extending ≳ 10 mm from the surface inlet) due to the slow chemical etching of the un-laser-irradiated material, but this can be compensated for by introducing inverse tapering in the specified design.
Long-rectangular channels can, therefore, be fabricated fairly easily; however, due to a slightly ellipsoidal laser spot (minor axis ≈ 2.6 μ m, major axis ≈ 6 μ m, Meineke et al. 2016), accurate circular channels are most easily obtained in the plane normal to the laser axis. Fabricating in this direction limits the maximum length of axisymmetric channel that can be etched in a single piece to around 5 mm (i.e., the substrate thickness). Accordingly, we scale our device geometries by setting R in = 0.36 mm, such that 13R in < 5 mm, allowing the central converging/diverging section of each geometry to be fabricated monolithically. We also fabricate 5 mm-long inlet and outlet sections of radius R = R in , which are then bonded either side of the central piece using ultraviolet curing epoxy resin. The three pieces are carefully aligned using locating pins passing through three corners of each individual section. Photographs of each device thus fabricated are shown in Fig. 2, with the dimensions labeled. Each photographed device is shown with the targeted wall profile superimposed in red, showing good fidelity in each case. The joins between the three individual glass pieces are visible in Fig. 2c, close to the dashed vertical lines that mark the limits of L opt . Note that to obtain these photographs, the channels were filled with a mixture of 82 wt.% glycerol and 18 wt.% water, which has a refractive index n ≈ 1.45 , close to that of the fused-silica substrate.
The fidelity of the fabrication to the specified design is especially critical for the optimized shape microchannel. To obtain a quantitative measure of the deviation between the fabricated channel and the specified shape, we performed an X-ray microtomography ( -CT) scan of the channel using a Zeiss Xradia 510 Versa 3D X-ray microscope operated with the Zeiss Scout-and-Scan Control System software. For this scan, the channel was filled with a water-based contrast agent (Omnipaque 350, Daiichi Sankyo Co, Japan, at a concentration of 350 g L −1 ), which provides a strong contrast against the fused-silica microchannel substrate. The 1601 projection scan was performed with a 4× magnification objective, 100 kV voltage, 9 W power, and 4 s exposure time. The -CT scan data were reconstructed using Amira analysis software (Thermo Fisher) and exported to Rhinoceros 3D modeling software (Robert McNeel and Associates) for rendering and comparison to the target profile (as presented in Fig. 3). Figure 3a shows a portion of the scan data rendered by Rhinoceros. The data allow the root-meansquare surface roughness of the channel to be estimated at RMS ≈ 0.7 μ m. Line profiles along the surface of the scanned channel taken at four equally spaced azimuthal angles are plotted in Fig. 3b in comparison to the target optimized profile. In general, there is clearly a very good agreement between the channel radius measured from the -CT scan and the specified channel design. There is an increased deviation from the design shape for |z| ≳ 2.3 mm, which marks the location of the joins between the three assembled glass sections. From measurements of the radius taken along the imaged channel surface at eight azimuthal angles equally spaced at = ∕4 , we estimate the average channel radius , and, thus, the deviation from the target profile TAR (z) = |R(z) − R TAR (z)| . This is presented in terms of a percentage of the local target radius R TAR (z) in Fig. 3c, (black symbols). Over the entire section of scanned channel, the deviation is ≲ 5% of R TAR . The maximum percentage deviation from the target profile occurs along the converging section of channel, close to the contraction throat ( z ≈ −0.6 mm), which equates to an absolute discrepancy of ≈ 5 μ m. The maximum absolute discrepancy from the target occurs in the outlet region z ≳ 2.3 mm, where the 3.5% error translates to around 12 μ m. We can also estimate the rootmean-square deviation of the channel from axisymmetry by comparing the measured radius at different azimuthal angles t o t h e l o c a l a v e r a g e r a d i u s . This is also shown as a percentage of the target radius in Fig. 3c (red symbols). Over most of the scanned length of the channel, the deviation from axisymmetry is < 1% . There are larger deviations in the inlet ( ≈ 1.7% ) and the outlet channels ( ≈ 3.5% ), indicating that the three individual glass parts were not precisely aligned during assembly, resulting in a small eccentricity. We are currently working to improve the accuracy of circular holes fabricated in planes parallel to the light path of the LightFab laser scanner. This will enable the production of longer axisymmetric channels in a single piece of glass, and will also allow the fabrication of more elaborate geometries with interconnecting axisymmetric channels.

Flow control
To create a flow loop, tubing connectors are bonded with two-part epoxy to the inlets and outlets of the glass microchannels and connected to 5 mL Gastight syringes (Hamilton, NV) via flexible silicone tubing. Flow is driven at a precisely controlled volumetric flow rate Q using two neMESYS low-pressure syringe pumps with 29:1 gear ratios (Cetoni GmbH, Germany), with one pump injecting fluid at the inlet and the second withdrawing from the outlet. The Newtonian fluid employed in the experiments is a mixture of 82 wt.% glycerol and 18 wt.% water, with a density ≈ 1210 kg m −3 and a dynamic viscosity ≈ 56 mPa s at the temperature of the experiments ( 25 ± 1 • C). The Reynolds number is a function of the position in the channel: Re = 2 Q∕ R(z) , which is maximal at the narrowest point, where R = R c . For the range of flow rates applied in the experiments ( 5 ≤ Q ≤ 500 L min −1 ), the corresponding Fig. 2 Micrographs of the three axisymmetric contraction/expansion geometries. a Abrupt, b hyperbolic, and c optimized. In each case, the specified wall profile is superimposed in red to indicate the accuracy of the fabrication Reynolds number at the contraction throat varies in the range 0.01 ≲ Re max ≲ 1 , so inertia can effectively be ignored.

Micro-particle image velocimetry ( -PIV)
The flow of the Newtonian fluid through the axisymmetric microdevices is quantified experimentally using a volume illumination -PIV system (TSI Inc., MN). The test fluids are seeded with a low concentration ( c p ≈ 0.02 wt.%) of d p = 2 m-diameter fluorescent tracer particles (Life Technologies) with emission wavelength = 575 nm. The microdevice of interest is mounted on the stage of an inverted microscope (Nikon Eclipse Ti) with the z-axis parallel to the stage (i.e., horizontal), and the cross-section of the flow geometry is brought into focus with an M = 4 magnification, NA= 0.13 numerical aperture Nikon PlanFluor objective lens. With this setup, the depth of the measurement plane (which is normal to the direction of light propagation) over which microparticles contribute to the determination of the velocity field is m ≈ 180 μ m (Meinhart et al. 2000). Excitation with a dualpulsed Nd:YLF laser (wavelength of 527 nm, time separation between pulses t ) induces particle fluorescence, and a highspeed camera (Phantom MIRO, Vision Research) operating in frame-straddling mode captures pairs of particle images in synchrony with the laser pulses. At each flow rate examined, the time t is set, so that the average displacement of particles between the two images in each pair is ≈ 4 pixels. All the flows in the present experiments are steady in time, so to reduce noise at each flow rate, 50 image pairs are processed using an ensemble average cross-correlation PIV algorithm (TSI Insight 4G). A recursive Nyquist criterion is employed with a final interrogation area of 16 × 16 pixels to enhance the spatial resolution and obtain 2D velocity vectors v v v = (v r , v z ) spaced on a square grid of 16 × 16 m. Further image analysis, generation of contour plots, and streamline traces are performed using the software Tecplot Focus (Tecplot Inc., WA).
Since the channels are axisymmetric and have a variable radius along z, the measurement depth, m , affects the apparent (measured) flow field to an extent that depends on z. For instance, at the inlet, m ≈ 0.5R in but at the contraction throat ( z = 0 ), m ≈ 2R c . Given the Poiseuille velocity profile across the channel radius, it can be anticipated that averaging vertically over the measurement depth will have an increasingly significant impact on the apparent flow field as the channel radius gets smaller. The contribution of microparticles to the determination of the velocity field at the measurement plane is weighted according to the distance r from the measurement plane (only in the direction of light propagation) as defined by the following function (Olsen and Adrian 2000;Bourdon et al. 2006): is an analytical expression derived by Bourdon et al. (2006), as an improvement to that derived by Olsen and Adrian (2000), for the determination of the relative weights of particles in a volume-illuminated -PIV setup using a cross-correlation algorithm to compute velocity vectors Fig. 3 X-ray microtomography ( -CT) of the optimized axisymmetric contraction/expansion geometry. a A section of the 3D-rendered -CT scan for z > 0 (diverging part). b Line profiles along the surface of the scan at four different azimuthal angles compared with the target (optimized) profile. Solid symbols represent z < 0 (converging region) and open symbols represent z > 0 (the diverging region). c Deviation of the channel profile from the target (black symbols) and from axisymmetry (red symbols), determined from measurements at eight equally spaced azimuthal angles and presented as the percentage of the local target radius in interrogation areas. The weighting function W (Eq. 4) is shown in normalized form in Fig. 4 alongside analytical Poiseuille profiles for flow in circular channels of radius R in and R c (normalized by the average flow velocity in each case, U). It is clear that in the inlets and outlets of the experimental microfluidic devices, where R = R in , the weighting function decays while there is only a minor variation of the expected flow profile, i.e., the measurement depth will not be expected to severely affect the determination of the flow velocity at r = 0 . However, at the narrowest part of the channel, where R = R c , there is a large variation of the expected flow profile within the spread of the weighting function; even particles located at the channel walls will contribute to the determination of the velocity field at the channel crosssection. In practice, Eq. 4 will be used in the present work to apply a depth-averaging to the numerical data that closely matches the theoretical depth-averaging of -PIV results, thus enabling the comparison between the experimental and numerical velocity fields.

Numerical optimization
The core of the numerical routine used to optimize the shape of the axisymmetric contraction/expansion was outlined by Pimenta et al. (2018) for planar flow-focusing devices. The routine starts with the generation of the geometry and mesh of the device, which has a deformable wall whose shape is controlled by a set of design points. After the flow has been numerically simulated, a cost function is evaluated over the axis of each trial geometry (set of design points) and fed as input to a black-box derivative-free optimizer (NOMAD, Le Digabel 2011), which generates a new trial geometry (set of design points) that attempts to reduce the cost function. These operations are repeated in loop until the evolution of the cost function drops below a pre-defined threshold, or after exceeding a pre-defined number of iterations. In this work, the algorithm described by Pimenta et al. (2018) was adapted to the contraction/expansion geometry, which required changes in the cost function definition and in the parametrization of the deformable wall.
Geometry parametrization The parametrized axisymmetric contraction/expansion is depicted in Fig. 5. The geometry has radius R in in the zone of constant cross-section, outside the optimized region. Point 'F' is held fixed at r = R in and z = −6.5R in . The 13 design points controlling the shape of the deformable portion of the bounding wall range from P 0 to P 12 , and form a Catmull-Rom interpolating spline. The axial coordinate of these points is fixed, whereas the radial coordinate of each point corresponds to a design variable. The design points are uniformly distributed over the axial direction between z = −6.5R in (point F) and z = 0 (point P 12 ). Note that P 12 , the last design point at z = 0 , has a modifiable radial location. Therefore, the contraction ratio is not imposed geometrically, but as will be seen next, the cost function is able to indirectly impose a pre-defined value based on mass conservation. Note also that the geometry is made symmetric in relation to plane z = 0 . For creeping flow  The radius is constant ( R in ) outside of the optimized region. Point 'F' is fixed in space, whereas points P 0 to P 12 form a Catmull-Rom inter-polating spline that controls the shape of the wall. The geometry is symmetric in relation to plane z = 0 of a Newtonian fluid, this ensures a symmetric flow regarding that same plane.
Cost function The contraction/expansion geometry was optimized with the purpose of obtaining an homogeneous uniaxial extensional flow, characterized by: in cylindrical coordinates with the z-axis aligned with the main flow direction. The extensional strain rate ̇ is positive in the contraction region and negative in the expansion region. However, obtaining an optimized device obeying this condition in the entirety of its volume is a difficult task (which, for an axisymmetric device, would probably require an annular stream of lubricating fluid to avoid the no-slip condition around the channel surface, e.g., Shaw 1975;Everage and Ballman 1978;Wang and James 2011). Rather than seeking to achieve a completely homogeneous flow field, it has been found in the case of 2D geometries more prudent to restrict the target of the optimization to the centerline of the device (Zografos et al. 2016;Pimenta et al. 2018;Haward et al. 2012;Galindo-Rosales et al. 2014), and we follow the same approach here for the 3D optimization of the contraction/expansion geometry. The flow becomes naturally extension dominated in the vicinity of the centerline as a consequence of the optimization procedure, but the extension is expected to be heterogeneous (i.e., with a non-constant strain rate) approaching the wall. Nevertheless, this effect is expected to be weaker in an axisymmetric geometry than in an equivalent device with a rectangular cross-section optimized in a single direction. This is because the Poiseuille flow profile in a circular duct is self-similar regardless of the duct radius, but the velocity profile in a rectangular channel depends upon the aspect ratio (Zografos et al. 2016). Furthermore, we reiterate that even though both the planar and the axisymmetric versions of this geometry are optimized only along the centerline, the two, nevertheless, provide fundamentally different modes of extensional deformation along that line. By limiting the target of the optimization to the axis of the device (i.e., r = 0 ), we, therefore, seek a device having a constant extensional rate ( v z ∕ z =̇= const ) in the contraction zone (defined as being over the range −z 3 < z < 0 ), and an equal but opposite extensional rate ( v z ∕ z =̇= −const ) in the expansion zone (between 0 < z < z 3 ). The expected velocity and strain-rate profiles along the axis of such device are shown in Fig. 6. The velocity is symmetric and the strain rate is anti-symmetric in relation to z = 0 . These are not essential conditions for our purpose, but they are desirable and they simply arise due to the symmetry of the geometry about the same plane and the creeping flow condition.
Following the approach of Zografos et al. (2016), a parabolic velocity profile (in the region z 3 < |z| < z 1 , see Fig. 6) was introduced to connect the region of constant velocity with the extension-dominated region. This smooth transition of velocity is arguably easier to achieve numerically than an abrupt transition, leading to smaller oscillations of the velocity as a function of the streamwise position (Zografos et al. 2016). The length of the transition region z = z 1 − z 3 should be sufficient to suppress the velocity oscillations, but should also be as short as possible to maximize the region of constant velocity gradient. From a small number of preliminary optimizations with varying values of z , the value of z = R in was selected as a reasonable compromise in this regard. The target axial velocity profile depicted in Fig. 6 is analytically described by: 3 + 2az 1 z 3 , a n d z 2 = (z 1 + z 3 )∕2 . The average and axial velocity in the upstream and downstream non-optimized region are denoted as U in and v z,in (for laminar flow inside a pipe, v z,in = 2U in ), respectively, and the axial velocity at the throat of the contraction ( z = 0 ) is v z,c . Note that, due to axisymmetry, the radial and azimuthal velocity components vanish on the flow axis, i.e., v r,theo (z) = v ,theo (z) = 0 at r = 0. Fig. 6 Velocity and velocity gradient profiles along the axis of the device as targeted by the numerical optimization. The velocity is constant for |z| > z 1 , follows a second-order polynomial in a transition region ( z 3 < |z| < z 1 ), and evolves linearly in the extension-dominated region ( |z| < z 3 ). The velocity profile is symmetric about z = 0 The numerical optimization was carried out for z 1 = 5R in , z 3 = 4R in and v z,c ∕v z,in = 16 . From mass conservation, this ratio of velocities imposes a theoretical contraction ratio of C r = 16 . This is a theoretical value, because, in practice, the radius of point P 12 , at the throat of the contraction, is a design variable, which, however, tends to R in ∕4 as the cost function is minimized.
The cost function is defined as: for each trial geometry, i.e., for each vector {r i } i=0...12 , where v z,num is the velocity profile sampled over the axis of the trial geometry and N is the number of sampling points (typically, N = 100 ), which are kept fixed in space independently of the geometry. Flow solver The velocity field was obtained from the solution of the steady-state Navier-Stokes equations for a Newtonian fluid in the Stokesian limit ( Re → 0 ). The equations were solved with a second-order finite-volume method, using rheoTool Alves 2016, 2017). A mesh dependence study was carried out to ensure the space convergence of the results presented in this work.

Results
In this section, we present the results of experimental flow velocimetry and finite-volume numerical simulations for viscous Newtonian fluid flow through the three axisymmetric contraction/expansion geometries. The experiments are carried out at low Reynolds numbers ( Re max ≲ 1 ) and the simulations are performed in the Stokesian creeping flow limit ( Re → 0 ). In all the results presented subsequently, lengths are normalized by R in and velocities are normalized by U in . Normalized quantities are specified by a superscript '*'.

Abrupt contraction/expansion
Results obtained from the abrupt contraction/expansion geometry are presented in Fig. 7. In Fig. 7a, we show a side-to-side comparison between the experimental (top) and numerical (bottom) velocity fields for creeping Newtonian flow over the whole contraction/expansion region. Note that in the experiment (here shown for Q = 5 L min −1 , Re max ≈ 0.01 ), due to the available field of view at 4× magnification, three separate measurements are required to capture the displayed range of z * . The three obtained fields are subsequently assembled manually to obtain the top half of Fig. 7a. Note also the different range of the color scales (representing the normalized velocity magnitude |v v v * | ) between the experiment and the simulation. This is due to the apparently lower experimental velocity in the narrow constriction, which is a result of the increased importance of the measurement depth in the lower radius section of channel (as explained in Sect. 2.3). In Fig. 7b, we show the profiles of the streamwise flow velocity v * z as a function of the radial position across the channel upstream of the contraction ( z * ≈ −8 ) and downstream of the subsequent expansion ( z * ≈ 8 ). The experimental data from different imposed flow rates (all Re max ≲ 1 ) collapse well, as expected, and the profiles compare well with numerical profiles for fully developed Poiseuille flow in a circular pipe. Figure 7c shows the profiles of v * z as a function of r * across the narrow section of channel close to z = 0 . Since the narrow section of channel has a constant radius R c , for smoothing purposes the profiles shown in Fig. 7c are, in fact, obtained by averaging the velocity over a section of the channel between −2.5 ≲ z * ≲ 2.5 . Here, the experimental profiles for different imposed flow rates also collapse well, but they are significantly lower than the numerical prediction. However, we find that by applying the weighting function (Eq. 4, Fig. 4) to the numerical data, we obtain a weighted numerical profile (shown by the dotted red line) that closely matches the experimental measurement. This shows that the true flow velocity in the narrow section of the axisymmetric microchannel is, in fact, in good agreement with expectations.
In Fig. 7d, we examine the streamwise fluid velocity along the length of the flow axis (along r = 0 ). The experimental data at different flow rates collapse as expected and, upstream of the contraction plane, the flow velocity has a constant value v * z,in ≈ 2 . Approaching the contraction plane, the flow accelerates sharply to reach a constant value inside the constriction. This is predicted to be v * z,c = 16v * z,in = 32 , due to the contraction ratio C r = 16 . In the experimental measurement, v * z,c appears to be rather lower than expected, but this is accounted for by the measurement depth of the -PIV setup, as shown by the weighted numerical prediction (dotted red line). On passing the expansion plane, the flow decelerates into the outlet channel back to a velocity v * z,in ≈ 2 , forming a curve that is rather symmetric about z = 0.
The streamwise velocity gradient (or extensional strain rate v * z ∕ z * =̇ * ) along the centerline axis is shown in Fig. 7e. Here, to reduce scatter resulting from taking pointwise spatial derivatives of noisy data, the experimental curve (blue dashed line) is computed from the average of all the experimental data shown in Figure 7d. The extension rate profile in the abrupt axisymmetric contraction/expansion geometry is characterized by a positive spike of width z * ∼ 1 , reaching a peak of ̇ * ≈ 80 at the contraction plane, followed by an anti-symmetric spike at the downstream expansion plane. Clearly, there is no lengthscale over which the extensional rate can be meaningfully described as uniform or homogeneous, which is a necessary requirement of an extensional rheometer. However, abrupt changes in the cross-sectional area of axisymmetric ducts are commonplace in real flow loops used in, e.g., biomedical and industrial systems involving viscoelastic fluids. It will be interesting and instructive to study the non-linear dynamics of such flows (e.g., McKinley et al. 1991;Rothstein and McKinley 1999) in an axisymmetric microfluidic system using modern visualization and measurement techniques.

Hyperbolic contraction/expansion
The results of the experiments and numerical simulations in the axisymmetric hyperbolic contraction/expansion geometry are summarized in Fig. 8. A side-to-side comparison between the experimental (top) and numerical (bottom) creeping Newtonian flow fields is provided in Fig. 8a. As for the abrupt geometry described previously, due to the limited field of view in the -PIV setup, the experimental field presented in Fig. 8a is an assembly of three separate measurements made at different z regions. Also note the reduced range of the velocity magnitude color scale for the experiment, which is due to the effect of the -PIV measurement depth. Nevertheless, there is obviously a rather good agreement between the experiment and the simulation. The apparent good agreement is confirmed by taking profiles of the streamwise flow velocity ( v z (r) ) at locations across the channel inlet and outlet (at z * ≈ −8 and z * ≈ 8 , Fig. 8b) and across the contraction throat ( z = 0 , Fig. 8c). The experimental data at different flow rates collapse well at the low imposed Reynolds numbers ( Re max ≲ 1 ). Across the inlet and outlet, the -PIV measurement depth has an almost negligible effect on the velocity measurement. The experimental profiles agree very well with the numerical prediction (Fig. 8b), indicating that the flow is fully developed before the start of the converging region and redevelops fully in the outlet channel. Across the contraction throat (Fig. 8c), the increased significance of the measurement depth causes a reduction in the measured flow velocity below the numerical prediction. However, as was the case in the narrow section of the abrupt contraction/expansion geometry, weighting the numerical prediction using Eq. 4 results in a better agreement with the experiment.
In Fig. 8d, we present the streamwise fluid velocity v z (z) along the centerline axis. The experimental data at different flow rates collapse as expected and are mostly well described by the weighted numerical prediction. The flow velocity at the contraction throat ( z = 0 ) is a little higher than given by the weighted numerical prediction (also evident at r = 0 in Fig. 8c), which probably indicates a small error in the channel dimension at that point. It can be estimated that the higher experimental flow velocity at z = 0 would be accounted for if the fabricated geometry had a radius R c ≈ 87.5 m at the contraction throat, i.e., just ≈ 2.5 m smaller than intended. We observe good symmetry of v z (z) about the plane z = 0. Figure 8e shows the streamwise velocity gradient ( v * z ∕ z * =̇ * ) along the axis of the hyperbolic channel. The derivative computed from the average of all the experimental data (dashed blue line) is shown in comparison with the numerical prediction (black dots), the weighted numerical prediction (red dots), and the nominal value of ̇ expected in this geometry, given by Eq. 3 (gray solid line). Considering the inevitable noise in the experimental profile, it is obvious that the fabricated geometry performs well in comparison with the numerical predictions. Approaching the hyperbolic region from upstream, the data show that the velocity gradient begins to increase ≈ 1R in before the start of the converging section and that the nominal value of ̇ * nom = 6.67 is only reached after a distance of about 1R in within the converging part of the channel. In the downstream expansion region, an anti-symmetric behavior is observed. Despite the clear discrepancy from the ideal behavior near the start/end of the hyperbolic-shaped converging/diverging regions, an approximately uniform extensional rate is achieved by the geometry over a substantial length either side of the contraction throat. The numerical simulation indicates that, within the hyperbolic part of the microchannel, ̇ remains within 1% of the nominal value between 0.29 ≤ |z * | ≤ 3.01 and within 5% of the nominal value between 0.22 ≤ |z * | ≤ 3.75.

Optimized contraction/expansion
In this subsection, we discuss the results of experiments and numerical simulations performed in the optimized shape contraction/expansion geometry, which was numerically designed with the intention to improve the homogeneity of the extensional flow field generated along the axis of the hyperbolic converging/diverging channel. A comparison between the creeping Newtonian flow field obtained experimentally (top half) and numerically (bottom half) is shown in Fig. 9a. A very good qualitative agreement is observed, subject to the same caveats given previously for the abrupt and hyperbolic contraction/expansion geometries.
As before, we extract experimental velocity profiles of v z (r) across the inlet ( z * ≈ −8 ) and outlet ( z * ≈ 8 ) channels (Fig. 9b), which, by comparison with the numerical prediction, indicate that the flow is fully developed both upstream and downstream of the optimized section of channel. Across the contraction throat ( z = 0 ), the experimental v z (r) profile is in reasonable agreement with the weighted numerical prediction (Fig. 9c), as was also the case for the abrupt and hyperbolic channels discussed above.
The streamwise velocity profile ( v z (z) along the centerline) in the optimized geometry is shown in Fig. 9d. In general, the experimental data are well described by the weighted numerical prediction. The measured flow velocity is slightly higher than expected at z = 0 , but the difference is accounted for by the small ( ≈ 1% ) discrepancy from the target shape at this location (see Fig. 3c). We also note some skewness of the experimental data about z = 0 , with slightly higher velocities at z < 0 than at z > 0 . This can also be accounted for by the deviation of the actual channel shape from the design, which, as can be seen in Fig. 3b,c, is slightly narrower than the target at z < 0 , but is closer to the target at z > 0. Figure 9e shows the streamwise velocity gradient ( v * z ∕ z * =̇ * ) along the axis of the optimized channel. It is once again clear that, notwithstanding noise in the experimental data, the fabricated device performs well in comparison to the numerical simulation. Furthermore, the numerical simulation follows the target profile of the optimization very well. In this case, we find that the velocity gradient remains within 1% of the nominal value ( ̇ * nom = 6.67 ) between 0.57 ≤ |z * | ≤ 3.48 and within 5% of the nominal value between 0.17 ≤ |z * | ≤ 4.05 . These approximately homogeneous regions of extensional flow are 7.0% and 9.9% (respectively) larger than those found in the simpler hyperbolically profiled channel. The effect of the optimization on the resulting streamwise velocity and velocity gradient profiles compared with the hyperbolic case is subtle and is not very clear by comparing Figs. 8 and 9. Accordingly, we present a direct comparison between the velocity gradient profiles obtained from the hyperbolic and the optimized geometries in Fig. 10. Here, for clarity, we only show profiles resulting from the numerical simulations compared against the nominal profile described by Eq. 3 and the target profile of the numerical optimization. Also, since all the profiles are anti-symmetric about z = 0 , we only show the converging region of the geometry to obtain a more enlarged and detailed view of the salient features. It is evident that the optimization introduces some oscillations into the velocity gradient profile, particularly upstream of the converging region where the large bulges appear in the wall profile of the optimized geometry. The main effect of the bulges in the optimized geometry is to delay the increase in the velocity gradient of the flow until further downstream compared with the hyperbolic case. The increase in the velocity gradient then occurs over a shorter length in the optimized geometry, more efficiently approaching the nominal value inside the converging region ̇ * nom = 6.67 . Similar observations were made by Zografos et al. (2016) for flows in optimized planar contraction and expansion geometries.

Discussion
In this work, we have examined viscous fluid flows through several axisymmetric microfluidic geometries featuring contractions and expansions, which can be useful for generating extensional flows and for modeling aspects of various biological and industrial fluid processing applications. Selective laser-induced etching (SLE) was used to fabricate the geometries in fused-silica glass and micro-particle image velocimetry ( -PIV) was used to experimentally measure the velocity fields inside the channels under conditions of laminar flow (Reynolds number, Re ≲ 1 ). We began by fabricating and testing a canonical microfluidic abrupt axisymmetric contraction/expansion, which performed in accordance with the predictions of finite-volume numerical simulations. Such devices are important due to their ubiquity in real-life flow loops, but do not generate a homogeneous extensional flow field (which is a key requirement for applications such as extensional rheometry). Next, we examined a contraction/expansion geometry with hyperbolic converging/diverging regions, which can be expected to provide a nominally uniform extensional rate anticipated to be useful for extensional rheometry. The fabricated device was again shown to perform as expected on the basis of numerical simulation and did, indeed, produce a uniform elongation rate over a significant proportion of the hyperbolic profile. Finally, we attempted to improve the hyperbolic shape by carrying out a numerical optimization of the geometry. The resulting complex axisymmetric geometry was fabricated by SLE and tested experimentally, again with good performance. The numerically optimized shape provided some improvement on the homogeneity of the extensional flow field compared to the hyperbolic geometry. However, the O(10%) increase in length over which the extension rate remained approximately uniform perhaps does not represent a substantial improvement in performance. In planar microfluidic devices, the advantage of numerical optimization techniques seems to be more evident in multiplestream devices such as T-junctions and X-junctions (e.g., Pimenta et al. 2018;Zografos et al. 2019;Haward et al. 2012;Galindo-Rosales et al. 2014), for which it is inaccurate to compute an analytical shape based on simple methods. The same is likely to also be the case in axisymmetric devices. In many cases, axisymmetric microfluidic geometries can provide a closer representation of a process than the 2D planar (rectangular) microchannels often used, but they are challenging to fabricate and to make measurements inside of. We have shown that SLE is a viable method for the accurate production of axisymmetric microdevices with complex shape profiles and we have demonstrated that optical probes (here -PIV) can be used to precisely measure the flow inside such devices. In this work, we have focused on developing axisymmetric microdevices designed to impose uniaxial and biaxial extensional kinematics to fluid elements. The radial compression and stretching of fluid elements differs fundamentally from the planar extensional kinematics generated in a rectangular channel. Our new devices may find uses in studying the deformation of fibers, drops, bubble, or cells in uniaxial extensional and compressional velocity gradients, where morphological dynamics will certainly differ from those observed in planar flows. The hyperbolic and optimized shaped axisymmetric geometries may also find applications in the non-linear dynamics and extensional rheometry of complex fluids, which are also expected to vary depending on the mode of extensional deformation. In future, we plan to study complex fluid flows (polymer solutions or suspensions) in these devices. We also intend to try optimizing other 3D geometries such as flow-focusing microdevices based on coannular channels, and ultimately Fig. 10 Numerical extensional rate profiles along r = 0 for the hyperbolic and optimized contraction/expansion geometries. The simulation results are compared against the nominal extension rate profile given by Eq. 3 and the target profile of the optimization procedure given by Eq. 6 the 6-arm 3D cross-slot device for generating uniaxial and biaxial extensional flows with a stagnation point (Afonso et al. 2010;Haward et al. 2019).