Probing flow-induced nanostructure of complex fluids in arbitrary 2D flows using a fluidic four-roll mill (FFoRM)

Engineering flow processes to direct the microscopic structure of soft materials represents a growing area of materials research. In situ small-angle neutron scattering under flow (flow-SANS) is an attractive probe of fluid microstructure under simulated processing conditions, but current capabilities require many different sample environments to fully interrogate the deformations a fluid experiences in a realistic processing flow. Inspired by recent advances in microfluidics, we present a fluidic four-roll mill (FFoRM) capable of producing tunable 2D flow fields for in situ SANS measurements, that is intended to allow characterization of complex fluid nanostructure under arbitrary complex flows within a single sample environment. Computational fluid dynamics simulations are used to design a FFoRM that produces spatially homogeneous and sufficiently strong deformation fields. Particle tracking velocimetry experiments are then used to characterize the flows produced in the FFoRM for several classes of non-Newtonian fluids. Finally, a putative FFoRM-SANS workflow is demonstrated and validated through the characterization of flow-induced orientation in a semi-dilute cellulose nanocrystal dispersion under a range of 2D deformations. These novel experiments confirm that, for steady state straining flows at moderate strain rates, the nanocrystals orient along the principal strain-rate axis, in agreement with theories for rigid, rod-like Brownian particles in a homogeneous flow.


A Fluidic Four-Roll Mill (FFoRM) For SANS Measurements
Without loss of generality, one can describe any 2D linear flow field by decomposing the local velocity gradient tensor into two parameters related to the deformation type and magnitude of the deformation rate at any material point in the flow 22 . Here E and Ω are the strain-rate and vorticity tensors. The flow type parameter is usually defined as (2) where |E| is the magnitude of the rate of strain tensor and |Ω| is that of the vorticity tensor (where = A A A ( : ) 1/2 ). For the approximately 2D, planar flows considered in this work, the flow type parameter takes values from −1 (pure rotation) to 1 (pure elongation) and describes the relative magnitudes of vorticity and strain-rate present in the flow. The strength of the flow can be described through the magnitude of the in-plane components of the velocity gradient tensor (Γ  ) defined as Γ = ∇  u . Hence, in 2D homogeneous flows the velocity gradient tensor can be expressed in terms of Γ  and Λ as [23][24][25] ScIENtIfIc RepoRts | (2018) 8 with axes defined in the velocity (1) and velocity gradient (2) directions, hereafter referred to as the in-plane axes. Written in this form, Γ Λ +  1 2 is equivalent to the shear rate (γ) for a simple shear flow (Λ = 0) and the extension rate (ε ) for the purely extensional planar flow (Λ = 1). The assumption that the flow is 2D, as assumed in equation (3), is approximately true in the microfluidic four-roll mill provided that the dimensions of the flow device in the out-of-plane direction are large enough that ▽u is dominated by in-plane contributions for the majority of the flow domain. Specifically, we shall see that the 2D approximation, namely that u 3  , is reasonable in the central region between the top and bottom boundaries (u 3 = 0 at the central symmetry plane), and near the central stagnation point.
The original microfluidic four-roll mill is sketched in Fig. 1a. It is similar to a cross-slot flow device, but is composed of eight channels where the relative volumetric flow rates in diagonally opposing pairs of four channels control the flow type within the device, while the average of the flow rates in these four channels controls the flow strength. The remaining four channels are held at constant (usually ambient) pressure. Although Lee et al. showed that the device can produce a reasonable approximation of the full spectrum of homogeneous 2D flows for Newtonian fluids, it is not clear whether it will produce homogeneous 2D flows when the fluid is non-Newtonian and/or viscoelastic.
In this work (summarized in Fig. 2), we explore this question, as well as the applicability of the fluidic four-roll mill (FFoRM) for SANS measurements. In particular, we show that a modified version of the microfluidic geometry of Lee et al. can produce tunable 2D deformation fields for at least some non-Newtonian fluids, and can be integrated with SANS instrumentation for measurement of the in situ microstructural configuration projected  onto the velocity-velocity gradient plane. CFD simulations were used to optimize the device geometry to produce flows for Newtonian fluids that are homogeneous within the scattering volume for flow-SANS measurements, and strong enough to deform complex fluid microstructures. After experimentally confirming the predicted flow behavior for a Newtonian fluid (glycerol), we then investigate the ability to produce nearly 2D flows with controllably variable flow types and deformation rates over a range of fluids with non-Newtonian rheological behavior using particle tracking velocimetry (PTV) measurements. Flow-SANS measurements are then demonstrated for arbitrary deformation types, and applied to measuring flow-induced changes in the microstructure of a dispersion of rod-like cellulose nanocrystals (CNCs) in deuterated water/glycerol (ratio 90:10).

Device Simulation and Parameterization
As indicated in the preceding section, we believe that the FFoRM geometry based upon the original design of Lee et al. has the best chance of satisfying the constraints for simple implementation and interpretation of flow-SANS measurements: namely a spectrum of homogeneous, 2D flows with sufficient flow strength and enough accumulated strain for the material to reach a steady state configuration.
To accommodate the flows of non-Newtonian fluids, we modified the Lee et al. design by removing the sharp protrusions from the channel dividers at the inlet to the center region. The elimination of these sharp protrusions and corners suppresses the tendency to produce non-axisymmetric or unstable flows due to the elastic instabilities associated with streamline curvature for viscoelastic fluids 26,27 . Furthermore, it also enables fabrication protocols that are less prone to machining errors that affect flow stability, while not compromising the device's ability to generate all flow types. These modifications allow for a simpler device geometry that can be parameterized by the scheme illustrated in Fig. 1. By specifying the flow rates in specific channels, the FFoRM can be controlled to produce arbitrary 2D flows for Newtonian fluids. This is schematically illustrated in Fig. 3 for the generation of extensional, simple shear, and rotational flows. Approximately, the average of the flow rates in the inlets and outlets (Q 2 + Q 1 )/2 determines the magnitude of the velocity gradient, while the ratio of these two flow rates (Q 1 /Q 2 ) determines the flow type.
Generally, devices with larger dimensions can generate flows that are spatially uniform over a larger area. However, they are limited in their capacity to generate flows at high strain-rates before fluid inertia influences the flow, as indicated by the Reynolds number (Re) which characterizes the relative magnitude of inertial contributions to the fluid momentum relative to the viscous stress Here, ρ is the fluid density, η its viscosity, and H is the inlet channel height. With respect to fluid flow, high aspect ratio (H/W) devices are ideal to eliminate 3D flow effects (e.g. large gradients through the thickness of the device from the upper and lower boundaries) when the possibility for creating slip surfaces (e.g. through a lubricating fluid layer or chemically/physically modified surfaces) is infeasible 28  . For the purpose of flow simulations, we only consider the impact of Re as this provides a physical constraint on the device operation. Additionally, the neutron beam spot is collimated to a 1 mm diameter circle that ensures adequate signal from the SANS measurement. In the final device configuration described later, we will show that the 3 mm device thickness ensures that the magnitude of in-plane gradients represents a large percentage of the total magnitude of gradients in the beam region of the device, thereby justifying the assumption of nearly 2D flow in the device. 2D CFD simulations were used to determine the impact of the geometric parameters shown in Fig. 1 on the deformations (rate and type) achievable and the uniformity of the flow in the center plane of the device near the stagnation point. Using COMSOL Multiphysics software, the Navier-Stokes and continuity equations (for a Newtonian fluid) were solved numerically for the velocity and stress fields in the device. In order to ensure the device design is sufficient for the previously stated objectives, we define our objective to minimize flow non-uniformity (which we define as the standard deviation of Λ within the beam diameter at the center-plane of the device) and maximize the deformation rate (defined as the average deformation rate in the center plane beam region) for a given value of Re at the flow types of interest. For the Newtonian fluid, we found that the maximum Re achievable before significant flow modification due to inertia was 0.2, Therefore, the impact of geometric changes were investigated at constant Re = 0.1 so that the deformation rates calculated are near the maximum achievable in the center plane of the geometry before inertial flow modification. In this work, we choose to limit our focus to the uniformity of the flow type rather than the uniformity of the deformation rate. Non-uniformities in flow type and deformation rate varied similarly with geometric changes, but with higher magnitude for variations in the flow type. With regard to non-uniformity in the flow type, the objective is to ensure that the device produces flows with a maximum standard deviation of center-plane flow type of σ Λ ≤ 0.1 for all Q 1 /Q 2 . Note that this is a lower bound estimate of the uniformity of deformation experienced by the fluid in the probing beam region, since as will be shown later the region of uniform flow may extend beyond the beam region in certain directions.
More detail on the device optimization can be found in the Supplementary Material Section B. After systematically varying all geometric parameters (C, W, D, and R) in increments of 0.5 mm, we found that the geometry that generates the highest deformation rates while maintaining σ Λ ≤ 0.1 for all Q 1 /Q 2 was a geometry with C = 5 mm W = 2 mm, D = 1 mm, and R = 1 mm. For the Newtonian fluid glycerol at 20 °C (with ρ = 1.261 g/mL and η = 1.412 Pa s), this device geometry generates maximum deformation rates (when Re = 0.1) of Γ ≥ −  s 10 1 for all flow types. For water at 20 °C (ρ = 1.00 g/mL and η = 10 −3 Pa s), the device generates maximum deformation rates Γ ≥ . −  0 01 s 1 for all flow types. For a circular region of homogeneous, linear flow, the average amount of strain accumulated at material points within the region can be numerically estimated assuming there is no strain prior to entering this region. The value varies depending on the flow type, ranging from 1.1 for purely extensional flows to 0.53 for nearly shear flows (Λ = 0.01), but is independent of the flow rate. Furthermore, the fraction of the beam that has accumulated a strain of 1 or more varies from 0.55 for an extensional flow to 0.83 for a near shear flow (see Supplementary Materials Section C). The strain accumulated in a flow with closed streamlines (Λ ≤ 0) is infinite. These estimations, based only on the flow within the beam, are conservative in the sense that any extension of the homogeneous flow outside the beam is neglected. The Lagrangian convection of microstructures from these regions into the beam region will contribute to the accumulated strain, and hence to greater uniformity of microstructural configurations within the beam compared to an estimate that assumes no strain is accumulated prior to entering the beam region.
In principle, a similar geometry design process could be used to optimize the device for other criteria (e.g. higher deformation rates or accumulated strain) or for other model fluids by utilizing different constitutive models in the generation of flow fields with CFD. However, in what follows we evaluate the device designed above for Newtonian fluids to assess its capabilities for generating stable flows for several types of complex fluids.

Experimental Evaluation of the Flow Fields For Newtonian and Non-Newtonian Fluids
The FFoRM device developed above is based on the ability of the flow device to generate homogeneous flows of Newtonian fluids within at least the scattering region centered about the central stagnation point. In this section, we use flow visualization experiments to characterize the flows that are actually generated in the device (Fig. 2b,c) for Newtonian and four ubiquitous classes of non-Newtonian fluids (shear thinning, yield stress, elastic (Boger), and viscoelastic). This investigation is intended to address utility of the current FFoRM geometry for several complex fluids, identifying those fluids for which the current design is satisfactory, and the challenges that are encountered for other types of fluids.
Newtonian fluid behavior. We begin with measurements for a Newtonian fluid. It is expected that the measured behavior will be accurately predicted by the numerical simulations from the design process, but comparison of measured and predicted data will provide some indication of our ability to experimentally determine the velocity gradient using particle tracking velocimetry (PTV). In addition, we utilize 3D simulations to evaluate any out-of-plane flows or out-of-plane velocity gradients that may develop for a Newtonian fluid as a result of the device's finite thickness (these will be assessed by experiments on the shear thinning fluid to follow). These 3D simulations were carried out using COMSOL Multiphysics software with similar specifications as the 2D simulations, but with added no slip boundary conditions on the upper and lower walls. The PTV measurements were made by seeding the test fluid with 10 μm hollow glass spheres in dilute (300 ppm) concentrations and capturing their in-plane trajectories in the center plane and at the stagnation point of the device. Using well developed algorithms for tracking tracer particles between frames and determining their velocity fields, we determine the spatial in-plane velocity fields and, subsequently, the in-plane velocity gradient fields by numerical differentiation with a locally weighted least squares regression 29,30 .
Streamlines at the center-plane of the flow device as predicted by simulation are compared to experimental streakline images generated by superposing PTV images from a single captured video. These results are included in Fig. 4 for three representative flows (extension, shear, and rotation) at constant Q 2 = 0.5 mL/min (Re = 0.0035). Likewise, simulated and experimental 2D flow type parameters (Λ 2D ) at the center-plane are spatially evaluated from the 2D velocity gradient tensor and are used as a scalar comparison of the full velocity gradient tensor.
White regions in the experimental indicate area where no particles were tracked, and noise in the experimental data can be attributed to tracking errors 29 as well as small (~4%) variations in fluid velocity through the finite focal depth of measurement (~0.3 mm). Experiments show no particles moving in or out of the depth of field, indicating that, near the center plane, fully 3D flows are minimal. We note that the location of the stagnation point can shift from the center when operating at conditions near shear flow (Λ 2D ~ 0). We attribute this shift to small variations in ambient pressure, which can drive a shift along the center line of the device (left to right in Fig. 4) where fluid velocity is near zero. Similar difficulties were reported in Lee et al. 21 and we further note that this shift in stagnation point does not noticeably change the flow type generated at the center of the device nor the long residence time in the beam region.
To quantitatively compare the flows of Newtonian fluids experimentally generated in the FFoRM with simulations, we calculate volume-averaged values of the flow type parameter and magnitude of the velocity gradient tensor for a region near the stagnation point in the center of the device. In the experiments, we choose a 1 mm diameter circular averaging region to match the region of uniform flow that the device was designed to achieve. In the simulations, the averaging volume is a 1 mm diameter, and 0.3 mm height cylindrical region around the nominal stagnation point in the center plane of the device chosen to match the focal depth of the experimental imaging system. Representative results are shown in Fig. 5 as a function of the ratio of volumetric flow rates Q 1 /Q 2 2 mm  for a constant value of Q 2 = 1.0 mL/min (Re = 0.007). As with the spatially mapped flow type parameter, the numerical values of Λ 2D and Γ  from simulations and experiments of a Newtonian fluid agree quantitatively, as we should expect.
Non-Newtonian flow behavior. Four representative complex fluids were formulated and rheological measurements under shear in a cone-plate rheometer were used to determine the steady state shear viscosity for a range of applied shear rates. Details of these formulations and measurements are included in the Methods section. Summarizing, the fluids were formulated with relaxation time scales that insured nonlinear rheological response for shear rates/ velocity gradients in the range less than 20 s −1 that can be realized in the flow device with viscosities in the range 0.1-5 Pa s.
Inelastic shear-thinning fluid. An aqueous dispersion of rigid, rod-like cellulose nanocrystals (CNCs) was used as a characteristic shear thinning fluid with minimal elasticity. CNCs derived from wood, such as those used in this study, have been found to be parallelepipeds with length, width, and height in the range of 100-200, 10-20 nm, and 2-5 nm respectively [31][32][33] . Such non-spherical particle dispersions exhibit flow-induced particle alignment in an imposed straining flow, which modifies their contribution to the fluid's viscosity (thinning in a shear flow and thickening in an extensional flow) 33,34 . Rheological measurements were used to identify 5 wt% CNC in 85.5 wt% D 2 O and 9.5 wt% glycerol as a model inelastic shear thinning fluid. The CNC concentration was chosen as the highest concentration before the emergence of nematic domains (indicated by the presence of birefringence at rest). The choice of solvent is a compromise between increased fluid-particle contrast in the SANS measurements and increased suspending medium viscosity (reducing the rotational diffusivity and thereby reducing the strain rates required for flow-induced orientation). The steady state shear response of the CNC dispersion is typical of a shear thinning fluid, and is characterized by a pseudo-Newtonian plateau in the viscosity at low shear rates with a zero-shear viscosity η 0 = 2.50 Pa s, followed by the onset of shear thinning at a shear rate of 0.5 s −1 (Fig. 6bi).
For such a semi-dilute suspension of orientable particles, the degree to which particles align is governed by the rotational Peclét number, where D r is the long-time collective rotational diffusivity of a particle in the dispersion, which can be related to that for a single, dilute rod 35 . Note that |E| is used as the characteristic deformation rate. Consequently, the alignment of the dispersion is dependent only on the rate of strain, and not the flow type. With these definitions, the onset of shear thinning occurs for Pe r = O(1). Thus, the shear rate corresponding to the onset of shear thinning (γ = 0.5 s −1 ) is used as an estimate the effective rotational diffusivity for this non-dilute system, i.e., D r = 0.5 s −1 .
Representative streakline images of flows, along with magnitudes of the experimentally observed flow type parameter and velocity gradient tensor, are reported in Fig. 6b for a nominal applied Peclét number of 4, where the shear thinning exponent is approximately n = 0.3. Under these conditions, we find qualitatively similar flow patterns to those observed for the Newtonian fluid (shown in Fig. 6a). The largest impact from shear thinning rheology is in the rotational flow case, where the deformation rate generated is reduced by half from the nominal (Newtonian) value. We hypothesize that this effect is due to thinning of the fluid in the entry channels which reduces the shear stress acting on (and tending to rotate) the fluid in the central region of the geometry. Other than this difference, the flows generated match the flows of a Newtonian fluid, suggesting that shear thinning alone has little effect on the ability to achieve arbitrarily variable flow types in the FFoRM at least for this degree of shear-thinning.
Yield stress fluid. A dilute, aqueous solution of Carbopol 934 was formulated for use as a model yield stress fluid. Aqueous dispersions of neutralized Carbopol 934 display yielding behavior 36 . From rheological measurements (Fig. 6ci), the value of the yield stress is quantified by the plateau in stress at low shear rates. The concentration of Carbopol 934 (0.2 wt%) was chosen so that the fluid exhibits a moderate yield stress (2.4 Pa) upon neutralization with 25 wt% KOH. It should be noted that the fluid is very strongly shear-thinning beyond the yield stress.
Streakline images of the carbopol solution in the FFoRM at low nominal deformation rates (Γ  nom ~ 0.05 s −1 ) are included in Fig. 6c. Qualitative differences from Newtonian behavior are noticeable in the streaklines for Q 1 /Q 2 values that produce shear and rotational flows for Newtonian fluids. Operation under nominal shear conditions (Λ 2D ~ 0) produces a flow that is more rotational (Λ 2D < 0). Streaklines for the rotational flow condition display a larger region of closed streaklines than for the Newtonian case. We attribute these differences in the flow fields to a more plug-like flow in the entry channels due to the strong shear thinning of the fluid near the confining walls. We find that flows can be tuned to produce the full range of flow types from extension to rotation, albeit with different operating conditions compared to a Newtonian fluid for the desired flow type that must be determined through velocimetry experiments (or, in principle, via simulations for a corresponding rheological constitutive model).
Purely elastic fluid. A polyethylene oxide (PEO) based Boger (purely elastic) fluid was formulated to study the impact of elasticity on flows in the FFoRM. Boger fluids are a class of elastic fluids that have a constant viscosity in shear flow over a range of deformation rates, but significant normal stress differences. Under steady flow conditions in a general flow, the strain rate relative to the fluid's elastic relaxation time (λ) is given by the Weissenberg number We formulated a water-based PEO Boger fluid such that it has a nearly constant shear viscosity (η = 0.08 Pa s) up to 100 s −1 and λ ~ 0.16 s (Fig. 6di) 37 . Streakline images of the elastic fluid flows for Wi = 2 are included in Fig. 6d and show significant flow modification from the nominal (Newtonian fluid) behavior. We find that, for flows nominally dominated by extensional deformation, this flow modification tends to generate flows with flow types closer to shear, presumably as a way to minimize extensional stresses. For flows nominally dominated by shear and/or rotation, flow modification begins upstream of the stagnation point as the fluid rounds the entry corner into the central region, which then significantly influences the behavior near the stagnation point. We note, however, that for nominally straining flows (Λ 2D > 0), the flow typically remains stable, i.e., it achieves a steady state that is invariant in time over the periods of measurement. By contrast, in the shear and rotational cases (Λ 2D < 0), the flows also tend to be unsteady, as indicated by streaklines that cross one another.
To further characterize the development of flow modification for the steady, nominally straining cases, we examine the development of flow modification with increasing nominal Wi for the purely elastic Boger fluid with Λ 2D ~ 1. Figure 7a includes streakline images from experiments at several values of the corresponding nominal Wi. Below Wi = 1, the flow is unmodified and matches that of a Newtonian fluid. As the nominal deformation  Figure. 7b includes the streamline images for the flow fields at the corresponding experimental conditions. We find that the Oldroyd-B model captures the flow modification observed experimentally. These results are consistent with earlier work, where it was reported that stretching type deformations of elastic fluids tend to be minimized due to large extensional stresses [40][41][42] . We thus conclude that the inability to generate extension dominated flows of elastic fluids when Wi > 1 is a physical limitation due to the fluid's rheology that remains as a significant challenge for the FFoRM type devices, or any device in which the fluid experiences an abrupt change in cross-sectional geometry.
Viscoelastic Fluids. A solution of CTAB/NaNO 3 wormlike micelles (WLMs) was formulated to investigate the impact of viscoelasticity (elasticity with shear thinning) on flows in the FFoRM 43 . Wormlike micelles are long, flexible chains formed by surfactants. Due to the noncovalent nature of the interactions that give rise to chains, in contrast to polymer solutions, wormlike micelles exhibit dynamic scission and reformation 44 . Despite this difference, the surfactant chains can entangle to form a viscoelastic network with a mechanical response that is similar to entangled polymer solutions. We formulated a 100 mM CTAB/300 mM NaNO 3 solution in water that shows a viscoelastic response (λ ~ 0.3 s) (Fig. 6ei) at ambient temperature 43 .
The observed flow behavior is shown for Wi = 2 in Fig. 6e. Qualitatively, the flows are similar to those for the elastic fluid. However, the tendency to produce a weaker (more rotational) flow is enhanced under conditions that would produce an extensional flow with Wi = 2.0 for a Newtonian fluid, with Λ 2D being only 0.36. Furthermore and similar to the purely elastic fluid, the shear and rotational flow fields are significantly corrupted and also appear to be unsteady.

FFoRM Operating Limits.
Given the flows produced for various rheological responses, it is important to summarize the limits of operation of the FFoRM device with respect to both maximum deformation rate and flow type placed by the fluid (including inertia or elasticity) as well as by equipment limitations, such as torque limitations of the pumps used to drive the flow (Fig. 8). Inertial limitations, which we find occur when Re > 0.2, can be overcome for a particular fluid by scaling down the size of the device. Torque limitations can be overcome by using pumps with a higher torque limit. These limits place upper bounds on the accessible space of Γ  and Λ 2D achievable in the FFoRM for any fluid (Fig. 8a). Furthermore, fluid elasticity leads to flows that are significantly modified away from the nominally applied flow set by Q 1 /Q 2 (when Wi > 1 and Re  1), or are unsteady in time (when Wi > 1 and Re > 1). These limits place further bounds on the accessible space of 〈Wi〉 and Λ 2D achievable in the FFoRM for fluids with appreciable elasticity (Fig. 8b). Efforts to overcome the inability to generate extension dominated flows of elastic and viscoelastic fluids will be the subject of future work, based upon different designs of the flow geometry.
We conclude that the FFoRM geometry in its current configuration has significant limitations when used with elastic or viscoelastic fluids. It can provide the full range of flow fields for inelastic shear thinning fluids and for yield stress fluids in the rheologically interesting range of deformation rates. This does not mean that it cannot be usefully applied for elastic and viscoelastic fluids, but it will be necessary to have detailed knowledge of the flow fields in order to interpret SANS-based data on the configuration of the microstructure for such fluids.

Flow-SANS Study of a CNC Dispersion Using the FFoRM
Following these initial investigations into the flow behavior of non-Newtonian fluids in the FFoRM, we seek to validate that the proposed FFoRM-SANS workflow (Fig. 2) yields interpretable results for a fluid whose microstructural response under flow is predictable. To accomplish this validation, we investigated the CNC dispersion that was previously described as a model inelastic shear thinning fluid, since the direction of orientation of nonspherical particles at low Pe r may be predicted even without the benefit of a rheological model. In this section, we discuss the full spectrum of experiments required for FFoRM-SANS to probe the orientation of the CNCs near the stagnation point in so-called strong flows (Λ 2D > 0), where we expect measurable microstructure orientation due to a combination of hydrodynamic and interparticle interactions 22 . In particular, we first measure the velocity and velocity gradient fields generated in the device using PTV for quantitative comparison to CFD simulations of a model inelastic shear thinning fluid (Fig. 2a-c). The measured microstructural response in the dispersion determined by SANS is then compared to asymptotic results for the orientation distribution of dilute rods in steady, homogeneous 2D flow to validate that the FFoRM achieves a measurable steady-state microstructural response (Fig. 2d).

CNC dispersion flow visualization.
In order to validate the ability to generate nearly homogeneous 2D flows of the CNC dispersion in the FFoRM, we now investigate in greater detail the fluid's flow uniformity and stability in the device. If uniform flows are generated in the probing beam region of the device, this provides confidence that FFoRM-SANS measurements will probe the CNC structure under the desired, uniform flow conditions, provided that the strain accumulated is sufficient for the dispersion to reach steady state. PTV was performed in a similar manner to the experiments for the Newtonian fluid, and velocity gradients were calculated for fluid flows near the mid-plane of the device. Furthermore, using this data, we seek to validate the use of a simple generalized Newtonian fluid (GNF) model -in this case, the Carreau model fit to the shear rheology data in Fig. 6bi -for numerical (CFD) simulations of the flow. The model, including best-fit values of the parameters, is reported in the method section. The CFD simulations were performed in COMSOL Multiphysics with the same 2D and 3D FFoRM geometries and corresponding boundary conditions. If successful, this validation indicates that CFD simulations can be used to predict operating conditions corresponding to a particular generated flow condition, circumventing the need for flow measurements at all possible flow conditions.
Representative results from the 2D CFD simulations and PTV experiments (Fig. 9) demonstrate that the FFoRM can generate extensional, shear, and rotational flows (and arbitrary combinations thereof) for the CNC dispersion well into the shear thinning regime. When evaluating the differences between the flows of a Newtonian fluid (Fig. 4) and the shear thinning CNC dispersion (Fig. 9), we find only minor differences in device operation (Fig. 10, solid and dotted lines). Notably, the shear thinning of the suspension decreases the maximum deformation rates achievable before Re = 0.2 by approximately half. As with the Newtonian fluid, we find excellent agreement in the flow type mapping between the Carreau fluid simulations and experiments with the CNC dispersion. This finding was not necessarily expected given that our choice of constitutive model does not show an extension thickening response that one would expect from microstructural models of rod-like suspension microstructure 45 . We note that this agreement is likely due to only moderate increases in the extensional viscosity as predicted for semi-dilute, hydrodynamically interacting rods of the dimensions and concentrations that we are investigating (~2 times that of a Carreau fluid for Pe r  1) 46 . As we did for the flows of the Newtonian fluid, we map the flow fields generated for the CNC dispersion under nonlinear flow conditions (Pe r > 1) in a 1 mm circular area near the stagnation point for quantitative comparison of generated velocity gradients in 3D CFD simulation and experiment (Fig. 10a). We find good agreement in both the generated deformation rate magnitude and flow type. As noted already, the addition of shear thinning only slightly modifies the operation of the FFoRM. Additionally, we measure the variation of velocity gradient through the thickness of the device under extensional flow conditions (Pe r = 3.5, Λ 2D = 0.93) to assess the 3D variation in the out-of-plane direction (Fig. 10b). We find that the simulations capture the variation in the magnitude of the strain rates through the thickness of the device. We note that the simulations predict a more parabolic variation in the strain rate through the device than the experiments, which indicate a plateau of nearly constant strain rate at the center of the device. We attribute these differences to particles being tracked from outside the focal plane of the measurement which would lead to a smoothing of the parabolic profile. From these measurements, we consider the simulations to provide a good representation of the flows generated for the CNC dispersion.   Furthermore, the experiments indicate that significant gradients along the device thickness (coincident with the SANS beam path) only occur in a region within 10% of the confining walls of the device. As such, the flow can be approximated as an effective 2D flow, which significantly simplifies the interpretation of the SANS data to follow. Moreover, with the fluid flow mapping of the CNC dispersion in the FFoRM geometry as well as a validation of the CFD simulations, we gain the ability of imposing arbitrary, uniform flows of the nanoparticle dispersion in the beam region of the device.

FFoRM-SANS measurements and data analysis.
We turn to validate the final step in the FFoRM-SANS workflow (Fig. 2d), measuring the microstructural response of the CNC dispersion using FFoRM-SANS. Validation of the FFoRM-SANS technique requires that we observe nanoparticle orientation under flow conditions consistent with previous studies of similar dispersions, including more recent experimental studies on CNC dispersions in shear flows 33,47 , and models of colloidal nanoparticle orientation under steady flows 45 . In this validation study, we seek to measure the absolute direction and strength (how probable the rods are to exist in the aligned conformation) of orientation as we vary the flow type and strain rate.
FFoRM-SANS measurements were performed to determine how the microstructure of the CNC dispersion responds under flow at moderate strain rates relative to their effective rotational diffusivity (as indicated by Pe r , equation (5)). Large changes in the orientation distribution of rod-like dispersions are expected in extension dominated ("strong") flows (Λ 2D > 0) in the nonlinear regime (Pe r > 1). To achieve these deformations in the FFoRM, we use the previously validated 3D simulations to calculate the operating parameters (Q 2 and Q 1 /Q 2 ) that yield a desired average range value of Pe r and Λ 2D through the beam region of the device. As mentioned previously, these values correspond to the average strain rates achieved in the scattering volume including the regions near the upper and lower walls of the device which will generate lower Pe r (<20% of the total scattering volume experiences less than half of the strain rate generated in the center plane, see Supplementary Materials Fig. S5), albeit with uniform Λ 2D . In this study, we performed SANS measurements where Λ 2D was varied from 0 to 1 and Pe r was varied from 0.02 to 50. Measurement of shear dominated flows at high Pe r was limited by inertial effects due to the higher volumetric flow rates required to generate shearing flows compared to extensional flows at similar Pe r .
For rod-like dispersions such as the CNC solution studied here, the orientational distribution of the particles is reflected in the degree and orientation of anisotropy in the two-dimensional scattering pattern (Fig. 11). In order to quantify the scattering anisotropy, we perform an angle-dependent annular average of the data in the q-range of 0.009-0.016 Å −1 to match previous rheo-SANS studies of similar dispersions 33 . The chosen range corresponds to a local maximum in the equilibrium scattering intensity, in which the scattering is dominated by interparticle correlations. As such, the degree of scattering anisotropy in this q-range reflects the degree of orientational order of the nanorod dispersion. To quantify this ordering, the angle-dependent scattering intensity was fit to a Maier-Saupe type distribution 48 , using nonlinear least-squares regression where φ is the azimuthal angle in the q x -q y (velocity/velocity gradient) plane with φ = 0 corresponding to the positive-q x axis and I max , α, and φ 0 are adjustable fitting parameters corresponding to the maximum scattering magnitude, amplitude of anisotropy, and the location of the minimum in the scattering intensity respectively, the latter of which represents the average orientation of the microstructure. We note that, as defined, the angle φ 0 describes the average orientation relative to the positive-q x axis in the detector plane, and is defined independently from any local coordinate frame of the flow in the FFoRM device. A representative scattering analysis is included in Fig. 11, and depicts the azimuthal averaging of 2D data and fitting to extract the orientation angle, φ 0 . For these anisotropic scattering patterns, we correlate the minimum in the scattering intensity to a maximum in the real space orientation probability distribution function, indicating the most probable particle orientation. Additionally, the alignment factor with respect to annular scattering intensity is calculated numerically using the trapezoidal rule via The alignment factor is commonly utilized in flow-SANS experiments to provide a measure of the degree to which the microstructure is aligned along the φ 0 direction 48,49 , and ranges from 0 to as high as 1 for isotropic and some perfectly oriented microstructures, respectively.
The modeling of dilute, non-spherical Brownian suspensions subject to arbitrary flows can be accomplished by solving a Fokker-Planck equation describing the orientation probability density function (OPDF) of particle conformations 50 . Microscopically, the effect of increasing strain rate is to preferentially orient particles along the principle strain-rate axis which competes with rotational diffusion due to Brownian motion. While an analytical solution only exists in the case of steady, irrotational flows, asymptotic approximations in the limits of small and large Pe r can provide insight into the expected behavior of the CNC in the intermediate Pe r regime 45 . In the limit of low Pe r , the first deviation in the OPDF from isotopic is for the rod-like particles to align along the principal strain-rate axis (see Supplementary Materials Section A). In the limit of large Pe r , particles align along the flow direction in planar flows other than pure shear. From streamline images of the flow, we can determine the absolute orientation of the outflow axis (φ outflow ). Then, assuming a region of uniform, linear flow within the beam (equation (3)), the absolute direction of the principal strain-rate axis (φ strain ) is given by This enables comparison of the orientation of the principle strain-rate axis φ strain in the scattering (q x − q y ) plane with the average particle orientation φ 0 determined from the anisotropic scattering patterns measured in FFoRM-SANS measurements on the CNC dispersion. Again, φ strain , φ outflow , and φ 0 are all defined relative to the q x direction to allow direct comparison for all applied flow types. Similar to mechanical four-roll mills, the absolute orientation of the principal strain-rate axis in the FFoRM does not vary appreciably with changes in the flow type parameter (Fig. 12), as confirmed through measurements of the outflow axis from simulated streamlines. Figure 12 shows representative results of this comparison for constant Pe r ~ 2. For reference, the outflow direction from the stagnation point measured with PTV is also shown. We find that for Pe r ~ 2, the CNCs align preferentially along the principal strain-rate axis for all flow types from 0 to 1, in agreement with the asymptotic prediction for low Pe r . This agreement between the orientation of CNCs and the principal strain-rate axis means that the current device configuration provides a sufficient amount of accumulated strain to orient the microstructure in the direction dictated by the local deformation field. The implications of this result will be discussed later. We further examine the dependence of the degree of alignment, as measured by A f , on various measures of the strength of the flow (Fig. 13). We find that A f increases monotonically with increasing deformation rate for conditions with similar flow type. However, this trend varies when plotted against 〈Γ  〉 (Fig. 13a). If instead the deformation rate is quantified in terms of Pe r = |E|/D r (i.e. the dimensionless rate of strain), the measured A f values exhibits modest collapse, with a trend that is independent of the flow type. This observation is consistent with the theoretical prediction for dilute, nonspherical particles at low Pe r, where the orientation distribution function (and therefore the alignment factor) is only a function of the Pe r . This result would be non-obvious in the absence of theoretical predictions, since the collapse of A f with strain rate is observed despite the larger upstream fluid velocities (proportional to Q 2 ) required to generate flows with lower flow type in the probing beam region.
From the FFoRM-SANS measurement of the CNC dispersion, we have extracted the direction of orientation and degree of alignment of a semi-dilute dispersion of CNC particles in the device. These results compare favorably to the asymptotic behavior of a model for dilute, rigid, rod-like dispersions in steady, 2D flows. The agreement between model and experiment suggests that, for the range of Pe r and Λ 2D measured, the orientation behavior of this semi-dilute CNC dispersion can be captured by defining the flow in terms of the dimensionless applied rate of strain, Pe r , which accounts for both the non-dilute concentration as well as the flow type. Additionally, this agreement suggests that the deformations applied in the beam region are sufficiently uniform such that the measured structure is similar to the expected steady state structure in a homogeneous flow. Thus, we conclude that the microstructural response being measured by SANS in the device represents a near microstructural steady-state, with little or no Lagrangian variation through the beam region. We attribute the achievement of near steady-state microstructural configuration on the fact that there is a stagnation point within the scattering region, where convection out of the region of uniform flow is slowed (i.e. the accumulated strain is high) compared to 'flow-through' geometries such as contractions and expansions.
To illustrate the importance of the stagnation point, Fig. 14 compares flow-SANS measurements when the probing beam is placed at different points in the FFoRM device, including locations near the stagnation point with large accumulated strains as well as those further from the stagnation point where material is rapidly convected into and out of the scattering volume. SANS spectra (Fig. 14a) correspond to measurements at various locations, flow types, and Pe r (Fig. 14b-d) within the FFoRM geometry. In this experiment we probed entry (1 and 9), exit (3 and 7), and transition (2, 4, 6, and 8) regions in addition to the previously reported stagnation point flow (5) under conditions corresponding to Pe r = 10 and Λ 2D = 0.95 in the beam region. From these results, the impact of the Lagrangian unsteady nature of the flow outside the stagnation point, as well as non-uniformity of flow type within these regions, is made apparent. Specifically, in several spectra (1, 4, and 7), the microstructure shows no preferred alignment (A f < 0.01) in any direction despite the presence of a sufficiently large local deformation rate (Pe r > 1) for extension-dominated (Λ 2D > 0) flows in the beam region. We attribute this observed lack of alignment to the lack of a preferred orientation direction in the average orientation distribution function when averaging CNC orientations from different flow histories. One spectra (3) shows a very high degree of alignment, despite the lack of strong flows in the beam region. We attribute the increased alignment to the convection of highly aligned microstructure from near the stagnation point along the outflow axis. Due to the high velocities in this region, the microstructure does not have enough time to adopt a near steady-state orientation distribution before being convected out of the beam region. This result clearly illustrates the significance of the FFoRM's ability to produce arbitrary flow types in the stagnation region of a nearly 2D flow for SANS measurements, as otherwise the interpretation of the resulting scattering anisotropy is convoluted by a number of effects including Lagrangian unsteady (i.e., history-dependent) responses of the fluid microstructure as well as large variations in flow type and deformation rate magnitude within the probing neutron beam.

Conclusions
In this work, we have reported the first flow-SANS sample environment and associated measurement methods capable of probing complex fluid microstructure under near-2D steady state flows with a range of variable flow types within a single device, the fluidic four-roll mill (FFoRM). Modification of a previously designed microfluidic four-roll mill geometry using 2D CFD simulations resulted in the ability to produce relatively uniform, steady state flows of varying flow type in fluids with both Newtonian and some non-Newtonian rheological responses. In cases where the flow can be easily controlled, the modified FFoRM enables unambiguous measurement and interpretation of microstructural changes within a probing neutron beam. Inelastic fluids (shear thinning and yield stresses) generate stable, Newtonian-like flows while fluid elasticity generally modified flows from their nominal behavior beyond Wi = 1, requiring further characterization of the flow fields of elastic and viscoelastic fluids for FFoRM-SANS studies. The resulting FFoRM-SANS workflow was validated through measurement of a shear thinning, rod-like Brownian dispersion. Quantitative agreement between PTV measurements and flow simulations verify the achievable flows within the device, and demonstrate that simulation can be used to predict the mapping between device control variables and the resulting flow. FFoRM-SANS experiments confirmed that the device achieves near steady state, spatially uniform microstructural responses when the probing beam is placed in the stagnation region of the flow. Furthermore, the measurements on rod-like suspensions confirm a number of theoretical predictions including the preferred particle orientation and scaling of the degree of interparticle alignment with both the type and strength of the applied deformation. These results provide a small glimpse of the new scientific insights and structure-property relationships enabled by the ability of FFoRM-SANS to make microstructural measurements of fluids under controllably variable complex flows.  Elastic (Boger) fluid. Polyethylene glycol (M v = 8000 g/mol, Sigma-Aldrich, CAS # 25322-68-3) as received was dissolved in filtered water at 37.5 wt% and mixed for one hour. Polyethylene oxide (M v = 4 × 10 6 g/mol, Sigma-Aldrich) as received was then dissolved in this solution at 0.08 wt% and mixed for two weeks with a magnetic stir bar at low speed (60 rpm) to avoid polymer degradation.
Rheology. The rheological tests were performed with an AR-G2 rheometer (TA Instruments) at constant temperature of 25 °C. To measure the rheological properties, steady-shear viscosity tests were performed using a Couette geometry for materials with low viscosity (cup diameter = 30.42 mm, truncation gap = 0.5 mm, and bob diameter and length of 27.95 mm and 42.15 mm, respectively) or a cone-plate geometry for materials with high viscosity and for normal stress measurements (diameter = 60 mm, angle = 2.00° and truncation gap = 0.055 mm). To prevent solvent evaporation, for the cone-plate, the material was thermally isolated with a solvent trap. Steady-shear viscosity tests were run in the range of shear rate 0.001 < γ <100 s −1 . Data were collected once every 300 s, in order to reach the steady state.
Computational fluid dynamics (CFD) simulations. CFD simulations were carried out with the COMSOL Multiphysics software package for inelastic fliuds or in OpenFOAM using the viscoelasticFluidFoam solver for elastic fluids 39 . Flow fields for the inelastic fluids were generated from simulations for a Newtonian fluid (ρ = 1.261 g/mL, η = 1.412 Pa s) and a generalized Newtonian fluid with viscosity described by the Carreau model where η 0 is the zero shear viscosity, λ is the relaxation time, and n is a parameter that determines the degree of shear thinning. Parameters were chosen as a fit of the steady shear rheology of the 5 wt% CNC dispersion, thus giving η 0 = 2.51 Pa s, λ = 2.06 s, and n = 0.3. A critical shear rate can be defined as the inverse of the relaxation time (γ critical = 0.49 s −1 ) and corresponds to the onset of shear thinning. The generalized Navier-Stokes and continuity equations were solved numerically for a meshed representation of the geometry using a finite element method. Constant volume flow rate boundary conditions were specified at the ends of the channels corresponding to Q 1 and Q 2 , constant pressure boundary conditions were specified at the ends of unlabeled channels, and no slip boundary conditions were specified for all other walls. Geometries were simulated in 3D when comparing to experimental data to ensure maximum accuracy and in 2D for geometry design to enable faster simulation times.
Fluidic four-roll mill (FFoRM) Fabrication. Devices consist of two outer plates squeezed around an inner plate containing the FFoRM geometry. The outer plates were constructed out of stainless steel and contain holes for the flow of temperature controlling fluid and holes for quartz windows to allow for light or neutrons to pass through the device. The inner plate was constructed out of titanium and the geometry was cut using a wire EDM. The central device geometry was cut according to the design in Fig. 1 and 35 mm channels were added as inlets/ outlets into the center geometry to allow time for the fluid to relax from entry effects. Holes were drilled though the sides of the center plate, connecting to the eight entry channels, to allow for fluid flow straight into the device.
Particle tracking velocimetry (PTV). Flows of the fluids in the FFoRM were experimentally determined using particle tracking velocimetry (PTV). The test fluid was seeded with glass spheres (10 μm diameter, 300 ppm) and injected into the device's four inlet channels using two syringe pumps (pump: Harvard Apparatus PHD 2000, tubing: Saint-Gobain Versilon 2001) at different rates. The four additional outlets were immersed in a container of the test fluid. The outlet tubing was cut to similar length and immersed at similar height to avoid differences in pressure between outlet channels. The device temperature was maintained at 25 °C by using a water bath connected to the device's outer plates. A CMOS camera (Point Gray Gazelle 2.2 MP Mono Camera Link) with 12x total magnification and variable frame rate (from 1 to 280fps) was focused on the center plane of the device geometry. The focal depth for the measurement is ~0.3 mm. Particle velocities were determined from videos using PTV algorithms that track individual particles with subpixel resolution and determine their velocities 29 . The velocity gradient tensor was calculated by fitting a second order surface using a weighted method of least squares 30,51 adapted for 2D. The search radius is an adjustable parameter and the fit is weighted using the tri-cube kernel: Flow-small-angle neutron scattering (flow-SANS). SANS measurements were performed using the NGB 10 m SANS instrument at the National Institute of Standards and Technology Center for Neutron Research (Gaithersburg, MD). A 1 mm diameter beam was collimated at the stagnation point of the geometry. Material of interest was flowed into the device using the same protocols as for PTV experiments. Measurements were delayed for 5 minutes after startup to ensure the flow had reached steady state. The scattering from the sample was collected in the q-range from 0.005-0.1 Å −1 with the wavelength λ = 6Å and a wavelength spread of Δλ/λ = 0.14.
The scattering vector q is defined as q = 4π sin(θ/2)/λ where θ is the angle at which the neutron is scattered and λ is the neutron wavelength. The temperature was maintained constant at 25 °C through the entire experiment with a water bath. The 2D scattering intensities in 128 by 128 channels were corrected for empty cell, plexiglass standard, and detector efficiency. Scattering spectra were reduced using standard NCNR protocols with Igor PRO software 52 .

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.