Extensional rheometry of mobile fluids. Part I: OUBER, an optimized uniaxial and biaxial extensional rheometer

Almost all flows of practical importance are comprised of both shearing and extensional kinematics. Prominent examples include flows through contractions or expansions, around obstacles and through branching junctions. Simple Newtonian fluids can be fully characterized by knowledge of their shear viscosity $η$ alone, since their extensional viscosity is known to also be constant and equal to $3η$⁠, $4η$⁠, or $6η$ for uniaxial, planar, or biaxial extension, respectively, where the coefficients 3, 4, and 6 are the respective Trouton ratio, $Tr$ [1,2]. By contrast, for viscoelastic fluids such as polymeric solutions and melts, which are widely present in industrial and biological processes, the situation is very different. Here, for sufficiently high strain rates, extensional flows are very effective at unraveling and orienting polymer chains [3–7]. Due to the entropic elasticity of the polymer, causing it to resist deformation, the hydrodynamically forced unraveling results in a non-linear increase in the elastic tensile stress difference and, hence, the extensional viscosity, with the nature of the increase being unknown a priori due to dependence on the fluid properties (e.g., polymer concentration, molecular weight, and extensibility) [8–11]. For viscoelastic fluid flows, even localized regions of extensional kinematics within the flow field can have a dominant impact on the macroscopic flow behavior [11]. For this reason, the quantitative characterization of the extensional viscosity of viscoelastic fluids is essential to enable a fully descriptive prediction of their behavior in arbitrary flow fields [12,13]. Unfortunately, the measurement of extensional viscosity is nontrivial, with a key challenge being to generate an extensional flow field that is both persistent and spatially uniform [2].

The potential of stagnation point flows for extensional rheometry has long been recognized, with the cross-slot device [Fig. 1(a)] being a relevant example. The cross-slot is a simple geometry, easily fabricated at microscale dimensions in order to obviate inertia with even low viscosity fluid samples. It consists of two oppositely facing rectangular inlet channels [aligned with $y$ in Fig. 1(a)] joined at right angles to two oppositely facing rectangular outlet channels (aligned with $x$⁠). By imposing an equal volumetric flow rate $Q$ through each of the four channels, an approximation to planar extensional flow is generated (ideally described by a rate-of-strain tensor $D$ with the only non-zero components being $D x x=− D y y= ε ˙$⁠, where $ε ˙$ is the strain rate). Crucially, there exists a free stagnation point at the center of symmetry of the cross-slot geometry where the flow velocity is zero. Hence, $ε ˙$ (⁠ $∝Q$⁠) is applied persistently, allowing strain to accumulate infinitely and any deformation of the microstructure (e.g., polymer) to reach a steady state. The steady state value of the tensile stress difference $Δσ( ε ˙)$ can be estimated by appropriate measurements of the pressure drop across an inlet and an outlet of the device [14], with the planar extensional viscosity given by $η P=Δσ/ ε ˙$⁠.

A problem with the “standard” cross-slot device as depicted in Fig. 1(a) is that the flow field is not homogeneous. The approximation to planar extension is only valid arbitrarily close to the stagnation point, and a given applied $Q$ does not provide a spatially uniform value of $ε ˙$⁠. By combining a finite-volume flow solver with an automatic mesh generator and an optimizer, Alves [15] iteratively modified the two-dimensional (2D) profile of the cross-slot in the region connecting the inlets and outlets in order to obtain an optimal approximation to ideal planar elongation [15,16]. The resulting optimized-shape cross-slot extensional rheometer [OSCER, Fig. 1(b)] imposes an almost homogeneous planar elongation over a region spanning $15W$ about the stagnation point, where $W$ is the characteristic channel half-width [15,16]. Since homogeneity is also required through the neutral $z$-direction, the flow should be 2D, so experimentally the OSCER geometry requires a high aspect ratio $H/W≈10$⁠, where $H$ is the channel half-height [see Fig. 1(b)] [16]. The OSCER geometry has proven useful for characterizing the extensional rheology and flow behavior of a variety of viscoelastic fluids [16–18].

Recently, advancements in three-dimensional (3D) microfabrication methods have motivated the development of a microfluidic six-arm cross-slot device [see Fig. 1(c)] [19,20]. Such a device can be operated in two modes. By injecting fluid at volumetric rate $Q/2$ along the two pairs of opposed inlets aligned with $x$ and $y$ and withdrawing fluid at volumetric rate $Q$ along the opposed outlets aligned with $z$⁠, an approximation to uniaxial extension (⁠ $D x x= D y y=− ε ˙/2$⁠, $D z z= ε ˙$⁠) is obtained. By reversing the flow in each channel, an approximation to biaxial extension (⁠ $D x x= D y y= ε ˙ B$⁠, $D z z=−2 ε ˙ B$⁠) is obtained. Note that the subscript “ $B$” on $ε ˙$ in the case of biaxial extension is to conform to established Society of Rheology notation [21–25] and distinguishes from the case where biaxial extension is considered as uniaxial compression, with $D x x= D y y= ε ˙/2$⁠, $D z z=− ε ˙$ [22,26]. The six-arm cross-slot device generates a stagnation point at its center and has been described as a microfluidic analog to the opposed jets apparatus [20,27,28]. Of benefit for extensional rheometry purposes, the microscale dimensions of the six-arm cross-slot minimize inertial effects in the flow, which were a severe limitation to the operability of the classical opposed jets apparatus [29,30]. Also, similarly to standard cross-slots, there is a possibility to evaluate $Δσ( ε ˙)$ (or $Δσ( ε ˙ B)$⁠) by appropriately designed pressure drop measurements. As such, the microfluidic six-arm cross-slot device possesses some promising attributes for use as a uniaxial and biaxial extensional rheometer for low viscosity “mobile” fluids. A question that arises is whether the 3D geometry of the six-arm cross-slot can be optimized in a way similar to the standard cross-slot device in order to obtain more homogeneous uniaxial and biaxial elongation.

Uniaxial and biaxial extension are kinematically the reverse of each other; uniaxial extension can equally be described as biaxial compression, while biaxial extension can be called uniaxial compression. Such flows have axisymmetry, with extension (compression) along an axis normal to the compressional (extensional) plane. With these considerations in mind, and given the practical constraints of a device with 4 (2) inlets and 2 (4) outlets, a first guess of the form of an optimized six-arm cross-slot is shown in Fig. 1(d). In such a device, to generate a uniaxial (biaxial) extensional flow, there are two opposing circular cross-section outlets (inlets) connected to four planar inlets (outlets). The inlets and outlets are connected by a nominally hyperbolic and axisymmetric shape, the precise optimal form of which must be determined.

In this work, we first perform a numerical shape optimization procedure on the six-arm cross-slot resulting in a number of 3D geometries that provide very close approximations to ideal uniaxial and biaxial extension under Newtonian creeping flow conditions. Selecting one of the optimal geometries for potential fabrication, we confirm its suitability for use with fluids of complex rheology by performing numerical simulations of the flow field with some commonly used viscoelastic fluid models. Next, we employ the 3D microfabrication technique of selective laser-induced etching (SLE) to fabricate a glass microfluidic geometry according to the selected design. Microtomographic particle image velocimetry (⁠ $μ$-TPIV) with a Newtonian fluid in the device, which we call the optimized uni- and biaxial extensional rheometer (OUBER), is used to demonstrate that the intended velocity fields are indeed imposed over a substantial volume of the device, generating a controlled and steady extension rate. Our ultimate aim with this research program is to undertake a comparison between the uniaxial, planar, and biaxial extensional rheology of viscoelastic fluids. It should be noted that although uniaxial and biaxial extension are the kinematic reverse of each other, the two flows are expected to result in significantly different polymer unraveling dynamics [2]. Accordingly, contrasting extensional rheology should result, although there is little convincing evidence from either experiment or theory to confirm or refute this speculation [23,25,31–33].

The remainder of the paper is organized as follows. In Sec. II, we describe the numerical optimization of the six-arm cross-slot geometry and we simulate viscoelastic flows in one of the resulting optimal geometries most suited to experimental fabrication. In Sec. III, we provide a detailed description of OUBER device fabrication and describe our flow measurement methods. The results of our Newtonian flow experiments in the OUBER device are presented and discussed in Sec. IV, and we draw our conclusions in Sec. V.

In this section, we will discuss the numerical optimization of the six-arm cross-slot device, originally proposed by Afonso et al. [19] and recently fabricated experimentally by Haward et al. [20]. To date, numerical optimizations of microfluidic device geometries have mostly been performed on 2D planar flow geometries, with the focus on achieving homogeneous extension rates [15,34–37]. Numerical optimization of a 3D flow geometry has only recently been demonstrated by Pimenta et al., who generated an optimized shape axisymmetric contraction-expansion geometry designed to produce a near constant extension rate along its axis [38].

In the present work, the numerical optimization of the six-arm cross-slot [Fig. 1(c)] follows a procedure similar to that outlined by Pimenta et al. [38]. Briefly, the geometry is parametrized with a set of design points that can be moved to deform the wall in order to achieve an extensional flow along a predefined region. The motion of these points is controlled by a derivative-free optimizer (Nomad v.3.9.12) [39], whose goal is to minimize an objective function embodying the difference between the flow in each candidate geometry and a theoretical homogeneous extensional flow. The velocity profiles used to compute the objective function are obtained with a finite-volume solver. These steps are discussed in more detail next and the interested reader can find more information about the automated optimization loop in Pimenta et al. [36,38].

The shape optimization of a six-arm cross-slot without any a priori consideration is a difficult task due to the large number of design parameters that can arise. However, as discussed in Sec. I, there is axial symmetry around the stretching direction in uniaxial extensional flows and around the compressive direction in biaxial extensional flows. Thus, in theory, a geometry with axial symmetry is required to impose such flows. This would naturally result in two opposing circular inlets (outlets) over the compressive (extensional) $z$-axis in biaxial (uniaxial) extensional flows that expand radially outward as they approach each other, to an outlet (inlet) over the $xy$ plane at $z=0$⁠. However, any experimental realization of such a device requires both inlet and outlet connections for the fluid flow. As apparent in our first guess of the shape of an optimized six-arm cross-slot geometry [Fig. 1(d)], the two opposing circular inlets (outlets) centered on the $z$-axis are not problematic; however, the four rectangular outlets (inlets) over the $xy$ plane inevitably break the axial symmetry of the geometry at some radial distance from the $z$-axis. This radial distance from the $z$-axis to the rectangular outlets (inlets) is parametrized in the optimization scheme as $L 1$⁠, while the distance along $z$ from the $z=0$ plane to the circular inlets (outlets) is $L 2$ (see Fig. 2). The radius of the circular inlets (outlets) is $R$⁠, while the half-width and half-height of the rectangular outlets (inlets) are $W$ and $H$⁠, respectively. This approach results in a geometry with several planes of symmetry, which can be obtained through successive reflections of 1/16 $th$ of the whole geometry [Fig. 2(c)], which, thus, represents the elementary unit to be optimized. Note that the unit element selected for optimization could be optimized on its entire 3D surface, without any assumption of axial symmetry. However, this option is not undertaken because a large number of design parameters would still be necessary to fine-tune the shape of the wall. Instead, we prefer to keep axial symmetry in the main body of the elementary unit and simply connect it to the rectangular cross section channel.

Under the approach described above, the wall of the elementary unit located between the circular and the rectangular cross section channels is parameterized through $n=13$ movable points, which form a Catmull–Rom interpolating spline [40], as depicted in Fig. 2(d). These points are evenly distributed through one quarter of an ellipse with center at $(x,z)=( L 1, L 2)$ and axes $( L 1−R)$ and $( L 2−H)$⁠. The radii of these design points (⁠ $r 1… r 13$⁠) can be adjusted by the optimizer to minimize the objective function, whereas their angle is kept fixed. The transition from the deformable section of the wall to the circular channels occurs at $z= L 2$⁠, whereas transition to planar channels occurs at $x= L 1$⁠.

Equation (1) corresponds to the formula adopted in this work to measure the objective function (⁠ $f o b j$⁠) for each candidate geometry $k$⁠, represented by its own array of design points ${ r 1 ⋯ n } k$⁠. Each of the two summations is relative to a given axis (⁠ $x$ and $z$⁠), as the velocity profiles are different along each one. For each candidate geometry $k$⁠, the velocity vector $u i , k$ (⁠ $u j , k$⁠), where $u=(u,v,w)$⁠, is sampled over $N$ (⁠ $M$⁠) points and compared with the theoretical or expected velocity vector $u i , t h e o$ (⁠ $u j , t h e o$⁠) at the given position $i$ (⁠ $j$⁠) of axis $x$ (⁠ $z$⁠). The constants $α$ and $β$ are used to weight the summation over each axis. Increasing the ratio of $α/β$ or $β/α$ increases the weighting of $f o b j$ over the $x$ or $z$ axis, respectively. In this work, we place equal importance on achieving homogeneous flows in both uniaxial and biaxial extension, so $α$ and $β$ are simply set equal to 1. Varying $α$ and $β$ in proportion has no effect on the final geometry but simply modifies $f o b j$ by the same factor. It should be noted that Eq. (1) corresponds to a simple, but effective formulation of a bi-objective optimization, as the motion of a given design point can improve the velocity profile over one axis but worsen the velocity profile over the other axis, i.e., there is a concurrent effect. More complex methods could be used in order to find the corresponding Pareto front (Nomad offers the possibility to do so) [39], but the simple single-objective formulation embodied by Eq. (1) is able to provide good results. This is significantly different from previous works [15,34–38], where a single velocity profile was imposed, either because a single axis existed [35–38] or because the velocity profiles over different axes were similar [15,34].

It can easily be shown that the average flow velocity across a circular cross section of increasing radius decreases with inverse proportionality to the radius. On the other hand, the rectangular channels impose a constant average and maximum velocity upon fixing $H$ and $W$⁠. Therefore, $L 1$ needs to be carefully chosen for each pair $(H,W)$ in order to avoid severe constrictions of the geometry over the $x$-axis. In practice, $(H,W)$ are selected first in the ratio range $3≤W/H≤4$⁠, and then $L 1$ is adjusted such that the velocity profile in a geometry with constant height $H$ over $x$ lies above $u x , t h e o$ (this indicates that $H$ should be increased in order to locally decrease the average velocity). Taking the velocity value at $L 1$ in such geometry, the value of $ε ˙/2$ is computed, which automatically defines the velocity profile over the $z$-axis, and hence the value of $L 2$⁠.

The governing equations are solved with the second-order finite-volume solver implemented in rheoTool [41],¹ which is based on OpenFOAM®. The geometry and mesh of the computational domain are built with standard tools provided in OpenFOAM®. A validation study was carried out to ensure mesh independency of the results obtained and presented in this work.

In this work, optimizations have been performed according to the scheme previously described using several combinations of the design parameters $L 1$⁠, $L 2$⁠, $H$⁠, and $W$ (see Table I), with each combination of parameters leading to a distinct geometry. In Fig. 3, the wall profile resulting from each imposed set of design parameters is shown, along with the corresponding axial velocity profiles obtained from the numerically solved flow field. In each case, there is an excellent agreement between the obtained velocity profile and the theoretical target.

The four optimizations of the six-arm cross-slot all yield geometries with similar performance, as evident from the similar respective values of the minimized objective function $f o b j , m i n$ given in Table I, and the close match between the numerical and theoretical axial velocity profiles in each case (Fig. 3). Hence, we only select one of them for fabrication and experimental verification. The most obvious suitable candidate geometry for the fabrication is geometry C, since it has the lowest value of $f o b j , m i n$ and it also does not possess the non-monotonicity in the wall profile that occurs in geometries A, B, and D at around $x/R=2$ (Fig. 3). Such deep concavities would not only be difficult to reproduce accurately by our fabrication method (see Sec. III A) but could easily trap air bubbles during fluid loading, affecting the resulting flow field during experimentation.

Since the ultimate intented application of the geometry is focused on extensional rheometry of viscoelastic fluids, it is important to assess the impact of fluid rheology on either the form of the geometry (or on the flow field imposed by the geometry). Note that, in principle, the geometry could be optimized using a viscoelastic rather than a Newtonian flow, as was performed in 2D for the planar OSCER geometry (in that case showing almost negligible differences in the resulting shape or velocity profiles) [15]. However, due to the large number of iterations (⁠ $∼$hundreds) needed to minimize the objective function [Eq. (1)], this approach is prohibitively computationally costly in 3D. Instead, we opt to demonstrate that our selected geometry (as optimized based on a Newtonian fluid) also imposes essentially the same flow field regardless of the rheology of the fluid. This is an important requirement to ensure that the device will be suitable for the characterization of different types of fluids. For this reason, we perform numerical simulations in geometry C using the Oldroyd-B and Phan-Thien and Tanner viscoelastic constitutive models, as described next.

The usual no-slip and no-penetration boundary conditions (i.e., $u= 0$⁠) are imposed on all surfaces of the channel. At channel inlets, we impose fully developed velocity and stress fields. At channel outflows, we apply the open boundary condition (OBC) [42]. Finally, we apply the usual symmetry conditions at all symmetry planes.

The Oldroyd-B (O-B) model is retrieved when setting $ε=0$ in Eq. (8), leading to $Y=1$⁠. Under steady simple shear, the O-B model predicts a constant viscosity, $η 0= η p+ η s$⁠, while the solvent-to-total viscosity ratio is defined as $β= η s/ η 0$⁠. Under steady extension (uniaxial, biaxial, or planar), the model predicts an extensional viscosity that is almost constant at low values of the strain rate, but which tends toward infinity as the dimensionless strain rate, or Weissenberg number, $Wi=λ ε ˙→0.5$ (⁠ $Wi=λ ε ˙ B→0.5$ in biaxial extension). We test two typically used values of $β$ (⁠ $1/9$ and 0.59) at two values of the Weissenberg number (⁠ $Wi=0.2$ and 0.4).

The linear version of the simplified PTT model (l-PTT) [43] is retrieved from Eq. (8) when $ε>0$⁠. The l-PTT model predicts shear-thinning effects in steady simple shear and a bounded extensional viscosity in steady extension (uniaxial, biaxial, or planar). With increasing $ε$⁠, the fluid becomes less strain-hardening and the onset of shear thinning is translated to lower values of the shear rate. We test typically used values of $β=1/9$ and $ε=0.02$ at two values of the Weissenberg number (⁠ $Wi=0.4$ and 0.8).

The Petrov–Galerkin stabilized finite element method for viscoelastic flows (PEGAFEM-V) [44,45] is used to solve the governing equations. We solve directly for the steady state solution, neglecting the time derivative in Eq. (6). The flow variables, $u$⁠, $p$⁠, and $τ$⁠, are interpolated by linear tetrahedra in a structured mesh. A validation study was again carried out to ensure mesh independency of the results at the values of $Wi$ examined herein (see  Appendix).

In Fig. 4, we present the results of viscoelastic flow simulations performed in geometry C. The upper left quadrant of Fig. 4(a) shows normalized velocity magnitude fields in the $y=0$ plane predicted for creeping Newtonian flow in uniaxial extension, while the remaining quandrants of the figure show the predictions of the viscoelastic fluid models at the highest values of the Weissenberg numbers tested. It is clear that the flow field predicted by the viscoelastic fluid models does not deviate significantly from the Newtonian prediction. Indeed, over the optimized region of the geometry, profiles of the streamwise velocity along the extensional axis (⁠ $w(z)$⁠) show excellent agreement with the Newtonian prediction for both the O-B and the l-PTT models at all values of $Wi$ [Fig. 4(b)]. Compared with the Newtonian velocity profile, we can observe only a slight overshoot in the velocity for the viscoelastic models near $z=±5R$⁠, and a slight reduction in the fully developed centerline flow velocity for the l-PTT model within the circular outlet channels for $| z |>5R$ (which is due to shear thinning).

Figure 4(c) shows normalized velocity magnitude fields in the $z=0$ plane predicted for flow in biaxial extension, again divided into quandrants depicting the creeping Newtonian flow (upper left) and the predictions of viscoelastic fluid models at the highest values of the Weissenberg numbers tested. Again there is visibly rather close agreement between the Newtonian and the viscoelastic predictions of the velocity field, although it is noticeable that the velocity magnitude tends to be slightly higher for the viscoelastic models close to the entrances of the planar outlet channels. Profiles of the streamwise flow velocity along the outlet axes in biaxial extension [i.e., $u(x)(≡v(y))$] again demonstrate rather close agreement between the predictions for Newtonian and viscoelastic flows [Fig. 4(d)]. The viscoelastic flow predictions deviate slightly from the Newtonian prediction near the limits of the optimized region (i.e., $x=y=±5R$⁠) and a small overshoot is evident immediately inside the planar outlet channels. Also, we can observe that for the l-PTT fluid, the fully developed centerline flow velocity within the outlet channels is slightly lower than for the Newtonian case (due to the shear thinning).

In general, numerical simulations indicate that even for viscoelastic models giving an unbounded response to extensional flow (i.e., O-B) and that account for the combination of elastic effects and a shear thinning viscosity (i.e., l-PTT), the flow field imposed by geometry C does not deviate significantly from that for the Newtonian creeping flow. The geometry provides the desired uniaxial and biaxial extensional flow fields given a variety of rheological conditions and imposed Weissenberg numbers up to $Wi=0.8$⁠. Therefore, it appears to be a promising candidate geometry for an extensional rheometer based on uniaxial and biaxial elongation.

In the remainder of the paper, we will focus on the experimental realization of geometry C and the verification of its performance based on the Newtonian fluid flow. The fabrication of the device, which from now on we will refer to as the OUBER, is achieved by the technique of selective laser-induced etching (SLE) in fused silica glass [46–48]. SLE is a two-step subtractive 3D printing technique for use with transparent substrates (typically glass). SLE in fused silica enables the fabrication of arbitrarily shaped microchannels with high resolution [ $∼ O(1μ$m)] in a rigid high modulus substrate with excellent optical clarity. In brief, the process involves the use of a scanning femtosecond laser to irradiate the volume to be removed (i.e., the internal volume of the microchannel) from a block of pristine fused silica substrate. The laser irradiation is performed using a commercially available LightFab 3D printer (LightFab GmbH). Subsequent to laser scanning, the fused silica block is ultrasonicated in potassium hydroxide at $80 °$C, and the irradiated material is selectively removed.

A 3D rendering of the design of the OUBER geometry, used to define the volume scanned by the femtosecond laser, is shown schematically in Fig. 5(a) (minus the inlet and outlet ports). The geometry is scaled such that the circular cross section channels have a radius $R=0.4$ mm. Thus, the half-width and half-height of the four planar channels are $W=0.64$ mm and $H=0.16$ mm, respectively (Table I). The geometry is divided into four parts along $z$⁠, each of which has a thickness $δz=5$ mm (the maximum thickness of substrate that can be used in the LightFab instrument). The reason for this division is because the laser scanning has a higher resolution in $x$ and $y$ than in $z$⁠. Therefore, circular holes are formed with the highest resolution in $xy$ planes. Scaling down the device dimensions to fit the entire channel within a single 5 mm thick substrate is not currently practical. The four individual pieces are assembled on locating pins and bonded together using ultraviolet-curing epoxy resin [see the exploded view of the assembly in Fig. 5(b)]. A photograph of the fully assembled glass device is provided in Fig. 5(c). x-ray microtomography (⁠ $μ$-CT) scanning of the central crossover region of the channel is performed using a Zeiss Xradia 510 Versa 3D x-ray microscope operated with the Zeiss Scout-and-Scan Control System software. The $μ$-CT scan data are reconstructed using Amira analysis software (Thermo Fisher) and exported to Rhinoceros 3D modeling software (Robert McNeel and Associates) to construct the 3D rendered image in Fig. 5(d). Eight surface profiles extracted from the $μ$-CT scan of the channel at four azimuthal angles (corresponding to positive and negative $x$⁠, $y$⁠, and $z$ directions) are compared to the target channel profile in Fig. 5(e). The inset to Fig. 5(e) shows the root-mean-square deviation $s z$ of the extracted profiles from the target (expressed as a percentage of the local target $z$⁠), demonstrating the excellent fidelity of the fabrication to the design (⁠ $s z≲0.05 z∀ x,y$⁠).

Note that the natural choice of the coordinate system has $x$ and $y$ axes aligned with adjacent planar channels and the $z$ axis aligned with the circular channels (as was done during the optimization step, cf. Fig. 2). However, since the axisymmetry of the device is broken by the four planar inlet/outlet channels, in the following, we will also consider a coordinate system rotated by $45 °$ about the $z$-axis [as shown in Fig. 5(a)], i.e., such that $x ′= 1 2(x+y)$⁠, $y ′= 1 2(y−x)$⁠. Comparison of the flow profiles along the standard and the $45 °$ rotated axes will reveal the extent to which the flow field remains axisymmetric.

Due to the surface curvature of the 3D OUBER device, see Fig. 5, clear imaging inside of the device (e.g., for performing flow velocimetry, as described below) requires that the channel be filled with a fluid of similar refractive index $RI$ as the fused silica glass. A sufficiently good match is achieved with a mixture of 89.6 wt. % glycerol and 10.4 wt. % water, with $RI=1.4582$ at $25 °$C (measured using an Anton-Paar Abbemat MW refractometer operating at 589 nm). This is close to the value of $RI=1.4584$ expected for fused silica under the same conditions [49]. The 89.6:10.4 wt. % glycerol:water mixture has density $ρ=1231 kg m − 3$ and viscosity $η=0.143$ Pa s.

Flow is driven through the microfluidic OUBER device by using 29:1 gear ratio neMESYS low pressure syringe pumps (Cetoni, GmbH) to control the volumetric flow rate through each individual inlet/outlet channel. For the uniaxial (biaxial) extensional flow, two pumps are used to impose a volumetric flow rate $Q$ through the two circular outlet (inlet) channels, while four pumps impose a volumetric flow rate $Q/2$ through the four inlet (outlet) channels. The pumps are fitted with Hamilton Gastight syringes of appropriate volumes so that the specified “pulsation free” dosing rate of each pump is always exceeded. Connections between the syringes and the microfluidic device are made using flexible Tygon tubing.

We consider the characteristic average flow velocity in the OUBER device as that in the circular cross section channels, $U=Q/π R 2$⁠. The Reynolds number of the flow is defined as $Re=2ρUR/η$⁠, and the maximum value reached in the experiments is $Re≈0.1$⁠, thus inertial effects in the flow are considered irrelevant. The expected extension rates obtained from numerical flow velocity profiles given in Fig. 3(c) are $ε ˙=0.4U/R$ and $ε ˙ B=0.2U/R$ in uniaxial and biaxial extension, respectively.

The flow field in the vicinity of the stagnation point of the OUBER device is measured volumetrically using microtomographic particle image velocimetry (⁠ $μ$-TPIV) [50]. Measurements are conducted using a LaVision FlowMaster system (LaVision GmbH), comprised of a stereomicroscope (SteREO V20, Zeiss AG, Germany) with dual high speed cameras (Phantom VEO 410, $1280×800$ pixels) imaging a fluid volume illuminated by a coaxial Nd:YLF laser (dual-pulsed, 527 nm wavelength). The fluid is seeded with $3.2μ$m diameter fluorescent particles (Fluoro-Max, Thermo Scientific), with excitation/emission wavelength of 542/612 nm, to a visual concentration of of $≈0.04$ particles-per-pixel.

Optical access to the stagnation point of the OUBER device is possible along the $y ′$ direction, between two of the planar channels (see Fig. 5). We focus on the $y ′=0$ plane of the OUBER device at $30×$ magnification, which enables reliable recording of the flow in a rectangular cuboidal volume defined by the limits $−1.4≤ x ′≤1.4 mm$⁠, $−0.6≤ y ′≤0.6 mm$⁠, and $−1.4≤z≤1.4 mm$⁠. The flow is recorded as single-frame images captured at a rate that is varied in proportion to the imposed flow rate such that no particle moves more than 8 pixels between consecutive frames. Images are preprocessed with local background subtraction and Gaussian smoothing at $3×3$ pixels. 3D calibration is performed by capturing reference images of a microgrid at the planes $y ′=±1000μ$m and $y ′=0μ$m, and a coordinate system is interpolated between these planes using a third-order polynomial. Particle positions in 3D are reconstructed from the images using four iterations of the Fast MART (Multiplicative Algebraic Reconstruction Technique) algorithm [51,52], followed by iterations of Sequential MART (SMART) [52], implemented in the commercial PIV software (DaVis 10.1.2, Lavision GmbH). We conclude the algorithm with five iterations of the sequential motion tracking enhancement (SMTE) method [53,54] to reduce the incidence of spurious “ghost” particles that arise due to randomly overlapping lines of sight [55] and which, thus, do not correlate in time. Volume self-calibration [56] is employed to improve the accuracy of reconstruction. Particle displacements between particle volumes are obtained using a multigrid iterative cross-correlation technique, with the final pass at $32×32×32$ voxels with 75% overlap yielding velocity vectors $u$ on a cubic grid of $33.2μ$m spacing. The obtained components of $u$ are labeled as $u ′$⁠, $v ′$⁠, and $w$ in $x ′$⁠, $y ′$⁠, and $z$ directions, respectively. Since the measured flows are time-steady, to reduce measurement noise, typically 50 vector fields are averaged (note that ghost particle intensity is converged after averaging of $≈5$ frames). In one particular case (for an imposed volumetric flow rate $Q=0.1$ mL  $min − 1$⁠), 200 vector fields are averaged in order to obtain sufficiently smooth data for computation of derived quantities. Subsequent to data acquisition, the software Tecplot 360 (Tecplot Inc., WA) is used for generation of contour plots, streamline traces, computation of vector components $u$ and $v$ (in respective $x$ and $y$ directions), for extraction of velocity profiles, etc.

In Fig. 6, we present experimental velocity magnitude fields (⁠ $| u |= u ′ 2 + v ′ 2 + w 2$⁠) with superimposed projected streamlines for a uniaxial extensional flow in the OUBER device at an imposed volumetric flow rate of $Q=0.1$ mL min $− 1$ (which corresponds to $Re≈0.02$⁠). The flow field as seen in the $y ′=0$ plane is shown in Fig. 6(a). The flow velocity decreases toward zero as the projected streamlines approach the stagnation point along the $x ′$ direction and increases with distance downstream of the stagnation point along the $z$ direction. Within the available field of view in the $y=0$ plane [Fig. 6(b)], the flow field appears similar to that in Fig. 6(a), as would be expected if the flow were ideally axisymmetric. In the $z=0$ plane [Fig. 6(c)], the field of view is quite restricted looking into the geometry along the $y ′$ direction; however, within the accessible field of view, we observe approximately circular contours of $| u |$ and streamlines that approach each other radially and converge at the stagnation point at $x ′= y ′=0$⁠.

Profiles of streamwise axial velocity components $u(x) | y = z = 0$ and $v(y) | x = z = 0$ are shown by open and closed symbols (respectively) for several imposed values of $Q$ in Fig. 7(a). Clearly (over the accessible field of view in $x$ and $y$⁠), the velocity profiles along the two orthogonal inlet axes are similar, with a nearly constant slope that becomes steeper as $Q$ is increased. Along the $x ′$ and $y ′$ axes [Fig. 7(b)], we have a different field of view in each direction. However, the axial profiles of the streamwise velocity components $u ′( x ′) | y ′ = z = 0$ (open symbols) and $v ′( y ′) | x ′ = z = 0$ (closed symbols) appear to be similar and to have an approximately constant slope as far as can be measured along $y ′$ (i.e., for $| x ′ |= | y ′ |≤0.6 mm$⁠). With increasing distance from the $z$-axis beyond $±0.6 mm$⁠, the profiles of $u ′( x ′)$ (open symbols) pass through local extrema before decreasing in magnitude. This is because the flow along the $x ′$ (and also the $y ′$⁠) axis is directed toward the boundary of the flow cell located at $L 1=5R(=2 mm)$⁠, where the flow velocity must vanish. Figure 7(c) shows profiles of the streamwise velocity component along the outlet axis $w(z) | x = y = 0$⁠. Over the measurable range of $z$⁠, the outlet axis velocity profiles are linear, with a slope that increases in proportion with $Q$⁠, as expected. Note that $u$⁠, $v$⁠, and $w$ all vanish at $(x,y,z)=(0,0,0)$⁠, which is the expected location of the stagnation point.

In Fig. 8, we present experimental velocity magnitude fields with superimposed projected streamlines for the biaxial extensional flow in the OUBER device, here again at an imposed volumetric flow rate of $Q=0.1$ mL  $min − 1$ (or $Re≈0.02$⁠). Figures 8(a), 8(b), and 8(c) illustrate the biaxial extensional flow field as observed in the $y ′=0$ plane, the $y=0$ plane, and the $z=0$ plane, respectively. Comparison with Fig. 6, for uniaxial extension at the same imposed $Q$⁠, shows that the velocity magnitude fields are almost identical, however the direction of the streamlines is reversed. This is as expected given the kinematic reversibility of uniaxial and biaxial extension. Accordingly, profiles of the streamwise axial velocity components in biaxial extension $u(x) | y = z = 0$ and $v(y) | x = z = 0$ [Fig. 9(a)], $u ′( x ′) | y ′ = z = 0$ and $v ′( y ′) | x ′ = z = 0$ [Fig. 9(b)], and $w(z) | x = y = 0$ [Fig. 9(c)] are essentially just mirror images of those obtained in uniaxial extension [Figs. 7(a)–7(c)].

Normalizing streamwise axial velocity components by the average flow velocity $U$ and normalizing distances by the radius $R$ of the circular cross section inlet/outlet channels, the experimentally measured axial velocity profiles for different imposed flow rates collapse, as expected for a Newtonian flow at low $Re$⁠. Mean normalized profiles computed from five measurements made for imposed volumetric flow rates $0.05≤Q≤0.8$ mL  $min − 1$ are shown for uniaxial and biaxial extension in Figs. 10(a) and 10(b), respectively. Note that these profiles are also computed by taking the mean of $u(x)$ and $v(y)$ and of $u ′( x ′)$ and $v ′( y ′)$ under the assumption that the flow along each of those two pairs of orthogonal directions is similar. The data points shown in Fig. 10 represent the normalized streamwise velocity profiles measured experimentally along the $x ′$ and $y ′$ axes (open circles), the $x$ and $y$ axes (closed up-triangles), and along the $z$ axis (open squares). The lines shown in Fig. 10 represent the target velocity profiles (i.e., the solutions of the Newtonian numerical simulations performed in the target flow geometry) along $x$ and $y$ (dashed line), along $x ′$ and $y ′$ (dotted line), and along $z$ (continuous line). Over the ranges of measurement, the experimental profiles clearly all agree very well with target numerical solutions. As mentioned above, if the flow were ideally axisymmetric, the profiles of $u(x)$⁠, $v(y)$⁠, $u ′( x ′)$ and $v ′( y ′)$ would all be identical. In Fig. 10 we observe that they agree well for $| x |≲R$ and $| x ′ |≲R$ (within $≈5%$⁠), but they diverge at greater radial distances from the $z$ axis, i.e., toward the perimeter of the circular region on the $z=0$ plane at $L 1=5R$⁠. For $| x |= | x ′ |=1.5R$⁠, $u ′≈0.85u$⁠, and for $| x |= | x ′ |=2R$⁠, $u ′≈0.75u$⁠. At $| x ′ |=5R$⁠, $u ′=0$⁠. Accordingly, in uniaxial extension [Fig. 10(a)], the measured extension rate along the $z$ axis is close to the numerical prediction $ε ˙=∂w/∂z=0.4U/R$ and is approximately uniform over the range $| z |≤5R$⁠. In biaxial extension [Fig. 10(b)], the extension rate over the $z=0$ plane is evidently also given by the numerical prediction $ε ˙ B=∂u/∂x(=∂v/∂y)=0.2U/R$ and is approximately uniform over a circular region defined by $x 2+ y 2≲ R 2$ (but is maintained over greater distances of $≈±5R$ along the $x$ and $y$ axes with which the planar outlet channels are aligned).

In Fig. 11, we compare the full field numerical prediction and experimental measurement of the extension rate in uniaxial elongation, $∂w/∂z$⁠, normalized by $U/R$⁠. Since the extension rate is a derived quantity, the experimental result is obtained from an average of mirrored and flipped velocity fields in order to smooth the data. Furthermore, we average data obtained from $x=0$ and $y=0$ planes and from $x ′=0$ and $y ′=0$ planes, on the assumption that the velocity fields over each of these two pairs of planes should be similar. In Fig. 11, the experimental result (only available in a limited field of view) is superimposed on numerical prediction over each imaged plane and is contained within the boundaries indicated by dotted gray lines. In all three planes { $x=0$ or $y=0$ [Fig. 11(a)], $x ′=0$ or $y ′=0$ [Fig. 11(b)], and $z=0$ [Fig. 11(c)]}, there is an excellent agreement between the numerical prediction and the experimental measurement, with closely matching contours of $∂w/∂z$⁠. Also, the normalized extension rate is close to the expected value of $ε ˙R/U=0.4$ over the large green regions observed in each plane.

In the same manner as Fig. 11, a full field comparison between numerically predicted and experimentally measured extension rates in biaxial elongation is presented in Fig. 12. Here, due to the two orthogonal axes of extension (along the $x$ and $y$ directions), we present the data in terms of $∂ u r/∂r$⁠, where $u r=ucos⁡θ+vsin⁡θ$⁠, $r= x 2 + y 2$⁠, and $θ= tan − 1⁡(y/x)$⁠. Once again, the extension rate is normalized by $U/R$⁠, revealing large green regions where the expected value of $ε ˙ BR/U=0.2$ is approximated and showing a generally good agreement between the experiment and the simulation in all three visualized planes { $x=0$ or $y=0$ [Fig. 12(a)], $x ′=0$ or $y ′=0$ [Fig. 12(b)], and $z=0$ [Fig. 12(c)]}. From these full field visualizations of $∂ u r/∂r$⁠, the difference between $x$ and $x ′$ (and between the $y$ and $y ′$⁠) directions is clearly evident, particularly from the view in the $z=0$ plane [Fig. 12(c)], where the expected strain rate is maintained for $≈5R$ along $x$ and $y$⁠, but becomes negative after $≈2R$ along $x ′$ and $y ′$ as the flow approaches the perimeter of the circular expansion region. The homogeneity of the flow field could very likely be improved by including four additional planar inlet/outlet channels aligned along the positive and negative $x ′$ and $y ′$ directions. However, this would add complexity to the experimental operation of the device, requiring additional syringes and pumps to control the flow and also further limiting the optical access for any desired quantification of the flow field.

Finally, to indicate the local flow kinematics, in Fig. 13 we present a comparison between the full field flow type parameter determined by simulation and experiment. The flow type parameter $ξ=( | γ ˙ |− | Ω |)/( | γ ˙ |+ | Ω |)$⁠, where $| γ ˙ |= γ ˙ : γ ˙ / 2$ is the magnitude of the deformation rate tensor, $γ ˙$⁠, and $| Ω |= Ω : Ω / 2$ is the magnitude of the vorticity tensor, $Ω$ [57]. Here, $ξ=−1$ indicates solid body rotation, $ξ=0$ indicates simple shear, and $ξ=1$ indicates purely extensional kinematics. Due to kinematic reversibility of uniaxial and biaxial flow configurations, the flow type parameter is expected to be the same for both (indeed the results obtained from the numerical simulations are identical). Therefore, in this case, we only show one set of fields and the experimental result shown in Fig. 13 is obtained by averaging the data from uniaxial and biaxial flows. Over the $x=0$ (or $y=0$⁠) plane [Fig. 13(a)] and over the $x ′=0$ (or $y ′=0$⁠) plane [Fig. 13(b)], there is an excellent agreement between the simulation and the experiment, with a very satisfactory matching between contours of $ξ$⁠. In the $z=0$ plane [Fig. 13(c)], the match between experiment and simulation is less impressive, with the experiment showing a reduced value of $ξ$ in comparison to the simulation. However, it must be remembered that the experimental result is derived from rather heavily smoothed and processed primary data. Also, velocimetry data on the $z=0$ plane has a somewhat low signal to noise ratio, lying directly along the line of sight into the flow cell and having a lower velocity magnitude than most of the field (see Figs. 6 and 8).

In general, it is evident that the regions of extensionally dominated flow kinematics (red regions where $ξ→1$ in Fig. 13) correspond to the regions of approximately uniform extension rates in uniaxial and biaxial flows (green regions where $∂w/∂z≈0.4$ in Fig. 11 and where $∂ u r/∂r≈0.2$ in Fig. 12, respectively). In short, the OUBER geometry successfully generates regions of nearly pure extensional flow at approximately uniform extension rate that extend over several characteristic device length scales in all three spatial directions. Taken as a whole, the $μ$-TPIV of the Newtonian flow field provides a clear confirmation that the fabricated OUBER device closely reproduces the numerically predicted flow fields and, therefore, has the potential for use as a uniaxial and biaxial extensional rheometer.

We have presented a numerical optimization of the “six-arm cross-slot” device [19,20] aimed at obtaining a geometry able to impose homogeneous uniaxial and biaxial stagnation point extensional flow fields with the intention of developing a uni- and biaxial extensional rheometer for mobile complex fluids. The optimization procedure (based on solving the Newtonian flow field) yielded a number of different geometries that depended on the input design parameters (i.e., the length scales over which the flow field was optimized). Of the generated geometries, one shape, in particular, was considered most amenable to fabrication and experimental verification of its performance. Prior to fabrication, it was confirmed by numerical simulations with Oldroyd-B and l-PTT models that the optimal flow field would also apply to rheologically complex constant viscosity and shear thinning viscoelastic fluids.

The device fabrication itself was achieved at microfluidic dimensions by the technique of selective laser-induced etching of fused silica glass. The fabrication yielded a highly precise match to the numerically designed geometry in a transparent substrate with optical access to the stagnation point region. Microtomographic particle image velocimetry for flow of a refractive index-matched Newtonian fluid at low Reynolds number was used to quantify the flow field in the experimental geometry over a relatively large volume centered on the stagnation point. These experiments provided confirmation of the good performance of the device.

In conclusion, we have designed, fabricated, and thoroughly tested a complex three-dimensional stagnation point microfluidic device, which has shown to produce good approximations to ideal uniaxial and ideal biaxial extension over multiple characteristic length scales in each spatial dimension. The applied extension rate scales linearly with the imposed volumetric flow rate, while the presence of the stagnation point means that the high (infinite) fluid strains requisite for steady-state extensional rheological measurements are achievable. Furthermore, the microfluidic dimensions of the device minimize the effects of inertia, which is essential for performing valid extensional viscosity measurements [30].

In Part II of this paper [58], we will demonstrate the use of pressure drop measurements in our new OUBER device for extracting the extensional rheological properties of viscoelastic fluids in uniaxial and biaxial extension. Furthermore, in combination with measurements made in the planar OSCER device [16], we will present a comparison between the uniaxial, planar, and biaxial extensional rheometry of model dilute polymeric solutions.

S.J.H., S.V., D.W.C., K.T.-P., and A.Q.S. gratefully acknowledge the support of the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan, along with funding from the Japan Society for the Promotion of Science (JSPS, Grant Nos. 21K14080, 21K03884, and 22K14184). F.P. and M.A.A. acknowledge the support provided by LA/P/0045/2020 (ALiCE) and UIDB/00532/2020 and UIDP/00532/2020 (CEFT), funded by national funds through FCT/MCTES (PIDDAC). We are indebted to Professor Robert J. Poole (University of Liverpool) for insightful discussions.

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.

In this appendix, we briefly outline the mesh dependence study carried out to ensure the results of viscoelastic flow simulations presented in Sec. II E are mesh independent. In all simulations, we use tetrahedral elements. To check the mesh convergence of our numerical solutions, we used three consecutively refined meshes, whose characteristics are quoted in Table II. Mesh M3 was used in all other simulations. The coarsest mesh (M1) is represented in Fig. 14.