Volumetric evolution of elastic turbulence in porous media

Abstract

Viscoelastic flow instability, which is compelled by elastic effects rather than inertia, can be driven to a chaotic state termed elastic turbulence (ET) manifested as strong velocity fluctuations with an algebraic decay in the frequency spectrum and increased mixing. We report the first spatiotemporally complete description of ET by considering a broad volume within a novel three-dimensional ordered porous medium, reconstructing flow at a micrometre characteristic length scale ($\text {Reynolds numbers} \ll 1$) via time-resolved microtomographic particle image velocimetry. Beyond a critical Weissenberg number of 2, we observe an elastic flow instability accompanied by an enhanced pressure drop with spectral characteristics typical of ET. Polymer chains in the ET flow state are advected along increasingly curved streamlines between pores such that they accumulate strain and generate a local flow instability evaluated per an established instability criterion based on local evaluation of elastic tensile stress and streamline curvature. The onset of ET leads to increased pore-scale resistance and positive feedback on upstream streamline curvature. ET is thus characterized by a continuous evolution between states of laminar and unstable flow: pores with unstable flow flood their adjacent peers and thus encourage straightened streamlines and flow stability across the array, while positive feedback from flow resistance on streamline curvature results in the instability propagating upstream along the array. By employing a geometrically ordered medium, we permit flow state communication between pores, yielding generalized insights highlighting the significance of spatial correlation and flow history, and thus provide new avenues for explaining the mechanisms of ET.

1. Introduction

Viscoelastic fluids exhibit both viscous and elastic behaviour owing to a microstructural component that responds elastically under strain, thus imparting an elastic response in addition to the viscous response. This class of fluids is commonly encountered in biological (e.g. blood and mucus) and industrial (e.g. polymer processing and surfactants) flows. Taking a polymer solution as a representative viscoelastic fluid, as the polymer chains advect over curved streamlines, they can accumulate stress while being elongated, and dissipate that stress over a characteristic time scale $\lambda$ while recoiling. A measure of the contribution of the elastic effect is the Weissenberg number ($Wi =\lambda \dot {\gamma }$, where $\dot {\gamma }$ is the shear rate). Beyond a critical $Wi_{c}$, accumulated elastic stress has been observed to drive purely elastic flow instabilities on curvilinear streamlines (McKinley, Pakdel & Öztekin 1996; Pakdel & McKinley 1996) despite negligible inertial effects; while at small Reynolds number $Re \ll 1$, Newtonian fluids would remain laminar and stable. The Pakdel–McKinley criterion has been used to capture effectively the onset of the elastic instability when the parameter $M = \sqrt {2\,Wi\,\lambda U/R}$ surpasses a critical $M_{c}$ dependent on the flow geometry and fluid rheology (McKinley et al. 1996; Pakdel & McKinley 1996). The $M$ parameter relates flow field perturbations and elastic stresses as a perturbation scale $\lambda U/R$ from the average flow velocity $U$ and streamline radius of curvature $R$, and takes $Wi$ as the strength of elastic tensile stresses. Thus a critical $M_{c}$ must be surpassed for an elastic instability to propagate in an arbitrary flow. The $M$ parameter has been shown to generalize the critical conditions for elastic stress buildup for a broad spectrum of flow geometries spanning microscopic to macroscopic scales, with a typical transition point of approximately $M_{c} = 1\unicode{x2013}7$ (McKinley et al. 1996; Morozov & van Saarloos 2007; Zilz et al. 2012; Haward, McKinley & Shen 2016; Qin et al. 2019; Kumar & Ardekani 2021).

At high $Wi \gg 1$ and low $Re \ll 1$, a chaotic elastic flow instability can arise bearing qualitative similarities to inertial turbulence (IT), thus termed elastic turbulence (ET). Both IT and ET exhibit substantial increases in flow resistance and mixing rate compared to their corresponding laminar regimes, and are comprised of a broad range of spatiotemporal modes with a spectral power-law decay exponent $\alpha$ (Groisman & Steinberg 2000, 2001; see review article Steinberg 2021). For ET, various studies have reported a steep $\alpha \approx 3.5$, with sensitivity to the specific flow geometry (Steinberg 2021). ET is thus temporally chaotic but spatially smooth, with dominance of a few low-order, vessel-scale modes, similar to IT flow on scales below the dissipation scale, viz. the Batchelor regime (Batchelor 1959; Burghelea, Segre & Steinberg 2005; Steinberg 2021). These characteristics lend utility to ET, and numerous industrial and environmental processes have been proposed that leverage the ET flow of polymer solutions through porous media for a spectrum of applications such as enhanced oil recovery (EOR), filtration, or groundwater remediation (Clarke et al. 2015; Howe, Clarke & Giernalczyk 2015). ET is driven by elastic stress reacting on the flow field, which is inherently difficult to quantify in situ. Furthermore, the spatiotemporal nature of ET means that a complete evaluation would be supported by the quantification of three-dimensional (3-D) flow fields spanning several vessel scales over long time scales. ET has been characterized in three dimensions primarily via direct numerical simulations for flows such as Taylor–Couette flow (Thomas, Sureshkumar & Khomami 2006; Thomas, Khomami & Sureshkumar 2009) and von Kármán swirling flow (Khambhampati & Handler 2020; van Buel & Stark 2022). Experiments have used holographic particle tracking to probe aspects of ET in three dimensions, such as particle dispersion (Afik & Steinberg 2017), while ET flow fields have been assessed about a 3-D cross-channel via particle tracking velocimetry and supported by pressure measurements (Qin et al. 2020). However, a volumetric quantification of ET near zero $Re$ evolving over a broad spatial scale (numerous vessel scales) remains lacking.

Regarding porous media flow in general, whereas viscous flow through porous media mixes poorly at the micro-scale due to the low $Re$, the elastic contribution provided by polymer stretching enhances the mixing in each pore and the overall dispersion of flow through the medium (Scholz et al. 2014). For shear-thinning solutions in this context, viscoelastic effects lead to a seemingly paradoxical pressure increase from the flow for increasing flow rates, despite the accompanying shear-thinning behaviour of the bulk solution (Marshall & Metzner 1967; Durst, Haas & Kaczmar 1981). Clarke et al. (2016) demonstrated that the viscoelastic effect in question was likely caused by elastic turbulence, supporting earlier observations in situ by Mitchell et al. (2016), and further validated the efficacy of ET for EOR. Recently, Browne & Datta (2021) were able to deduce an energy balance to relate directly the energy dissipation of a random porous medium to the pore-scale flow fluctuations of ET above $Wi_{c} = 2$. However, they noted that the energy balance neglected extensional viscosity and strain history effects, and that these effects likely became significant for $Wi > 4$ where the energy balance underpredicted the losses observed from experiments.

In a prior work, Browne, Shih & Datta (2020a) investigated viscoelastic flow through one-dimensional (1-D) model porous media for a system scale ranging from a single pore up to a thirty-pore sequence. Their experiments highlighted the importance of spatial correlation. Whereas flow above $Wi_{c}$ through a single pore generates an upstream vortex typical of a viscoelastic contraction flow instability (McKinley et al. 1991; Rothstein & McKinley 1999; Carlson, Shen & Haward 2021), the presence of sequential pores results in a reduced $Wi_{c}$ and the occurrence of bistable states of laminar and unstable flow. Flow through pores under a laminar state exhibits a corner vortex, while unstable flow is typified by strongly fluctuating streamlines flooding the corners of the pore. They hypothesized that the flow states were driven by competition between the extension and relaxation of polymer chains as they advected through the array. This notion was supported by later numerical work by Kumar et al. (2021). For closely placed pores in a 1-D array, Kumar et al. showed that the laminar flow through pores with eddy-dominated flow was accompanied by increased polymeric stress (extended polymer chains), whereas unstable flow through thus eddy-free pores exhibited reduced stress (coiled polymer chains).

Of relevance to these findings is the time scale ratio for which they manifest. Within the context of repeating pores, the pertinent time scales can be related as a pitch-based Deborah number, generally defined as $De = U\lambda / L$, where $L/U$ provides a time scale for material advection between pores separated by distance $L$. For $De > 1$, microstructures are advected between pores faster than they can relax. Kumar et al. (2021) noted that the flow instability for a ten-pore sequence started at $Wi > 18$ and $De > 1$: pore interaction affects the onset of instability for sufficiently large arrays. A later work reported a similar onset with $De$ for an instability occurring between two cylinders (Kumar & Ardekani 2021), and experiments by Ekanem et al. (2022) found that increased $De$ lowers $M_{c}$ for flow instability about separated pores. Indeed, the significance of $De$ for flow instability in repeating arrays was shown in experiments by De et al. (2017) for a two-dimensional (2-D) ordered array of cylinders. They found coexisting lanes of fast and slow flow in the array above $De = 1.5$, leading to flow crossing between channels driven by strong crossflow action. They conclude that $De$ is more apt than $Wi$ to indicate large-scale instabilities in repeating arrays. Walkama, Waisbord & Guasto (2020) and Haward, Hopkins & Shen (2021) considered viscoelastic flow through ordered and increasingly disordered cylinder arrays, and both works show that ET is likely controlled by the concentration of stagnation points in 2-D arrays. This is an analogous explanation for the $De$ scaling by De et al. (2017), where they in effect also modulated the concentration of stagnation points while geometrically scaling cylinder spacing for different $De$.

Experiments for viscoelastic flow through porous media encounter common limitations: the media are typically opaque, and flow is spatiotemporally sensitive to the specific geometric configuration of the medium (Walkama et al. 2020; Haward et al. 2021; Browne & Datta 2021; with reviews by Browne, Shih & Datta 2020b; Kumar, Guasto & Ardekani 2022). Refractive index matching (RIM) such that the medium and flow are transparent permits direct optical access for image-based interrogation techniques. Such methods include streak imagery and particle image velocimetry (PIV) for 2-D quantification, while recent works probed 3-D flows via particle tracking velocimetry (PTV) (Holzner et al. 2015; Souzy et al. 2020) and tomographic PIV (TPIV) (Larsson, Lundström & Lycksam 2018; Forslund et al. 2021). Souzy et al. (2020) used scanning PTV to resolve Newtonian creeping flow through randomly packed spheres (2 mm in diameter), achieving a fine spatial resolution for a time-averaged grid to show that 3-D mixing results in Lagrangian chaos for material advection. This result expanded upon prior numerical works (Turuban et al. 2018, 2019) that showed similarly enhanced mixing for flow past a body-centred cubic array of spheres compared to a simple cubic array.

For TPIV, synchronized cameras image a region of interest (ROI) from different lines of sight (LOS) and reconstruct the position of seeded particles within a volume of interest (VOI) for local cross-correlation in time to yield a velocity field. Applied for micro-scale, a stereo microscope provides limited optical LOS such that $\mathrm {\mu }$-TPIV can be conducted: recently, we have shown this approach to be useful for capturing transient viscoelastic flow instability in three dimensions (Carlson et al. 2021). While we do not match the spatial resolution of scanning PTV, $\mathrm {\mu }$-TPIV can time-resolve flow, which is crucial for quantifying the unsteady volumetric velocity fluctuations of ET.

In the present work, we apply $\mathrm {\mu }$-TPIV to the study of ET in a micro-scale ordered 3-D porous medium. Besides the applied virtues for studying viscoelastic flow through porous media, we can probe a variety of length scales by imaging simultaneously multiple pores and throats across a constant focal plane, thereby accessing the flow history between pores. Thus this configuration enables a more fundamental study of the spatiotemporal evolution of ET, where due to the dominance of vessel-scale (i.e. pore-scale) modes, ET flow will appear spatially smooth at a given pore for a temporal instant. The desired geometry is therefore a geometrically ordered fully 3-D array of cubicly packed glass spheres – a hitherto unstudied configuration, with fabrication made possible via selective laser etching (Gottmann, Hermans & Ortmann 2012; Meineke et al. 2016; Burshtein et al. 2019). We utilize $\mathrm {\mu }$-TPIV via RIM to assess flow through this novel 3-D ordered porous medium, together with pressure measurements, to provide a spatiotemporally complete quantification of ET.

From these experiments we scale the transition from laminar, to elastically unstable, to ET flow against the $M$ parameter, showing via localized ($M_{L}$) calculation over pore-scale instantaneous streamlines that coexisting regions of unstable and laminar flow hold locally consistent with the $M$ instability criterion. In the ET state, which we validated via comparison against common spectral characteristics, we relate the seemingly disparate coexisting states of laminar and unstable flow as a continuous progression between the two that comprises the core velocity fluctuations of ET. Finally, we relate inversely correlated flow states via competing streamline curvature, and the spatial evolution between states, as positive feedback between flow resistance and upstream streamline curvature.

2. Experimental methods

2.1. Flow cell and refractive index matching

The test cell was designed to confine flow about a $70 \times 5 \times 5$ array of $500\,{\mathrm {\mu }}{\rm m}$ diameter cubicly packed solid spheres (spaced $w_p = 500\,{\mathrm {\mu }}{\rm m}$) placed in a square-sectioned channel of width $w_c = 2.5\,\mathrm {mm}$, yielding nominal porosity $\phi = 0.4675$. The array nature of the geometry permits a simultaneous study of numerous pores at a given focal plane, and due to the constant $De$ throughout the array, we can assess the local flow history from time-resolved measurements. To fabricate this system, we used selective laser-induced etching (SLE) of an open-faced slab of fused silica (${\rm SiO}_2$) (Gottmann et al. 2012; Meineke et al. 2016; Burshtein et al. 2019) on a commercial LightFab SLE printer (LightFab GmbH, Germany), followed by chemical etching in potassium hydroxide to remove unwanted material. The open-faced flow cell was subsequently sealed with glass slides about the $z$ direction. As manufactured, the spheres are slightly over-etched and remain connected by glass ‘necks’ (images in figure 1a). To incorporate the actual sphere contacts in a stereolithography (STL) surface for $\mathrm {\mu }$-TPIV masking, we used X-ray microtomography ($\mathrm {\mu }$-CT; Pimenta et al. 2020) to quantify precisely the surface of the geometry (shown in figure 1b). Flow was driven through the device with a volumetric flow rate $Q$ via two syringe pumps (neMESYS, Cetoni GmbH, Germany) in a push–pull configuration, with a wet/wet differential pressure transducer (OMEGA, USA) spanning the inlet and outlet ($p_1$, $p_2$ in figure 1c) sampling at 100 Hz. The characteristic velocity of flow passing about the array of spheres is $U_{p} = Q/A\phi$, the average flow velocity in each pore between spheres for a volumetric rate $Q$, total cell cross-sectional area $A = w_{c}^2$, and porosity $\phi$.

Figure 1. (a) Images of the model porous medium as manufactured and filled with water for contrast.(b) A section of the micro-CT of the porous array together with the alignment of the volume of interest (VOI, blue), centred 35 pores into the array. The fluid–solid interface is coloured white. (c) The beam path of the LaVision FlowMaster $\mathrm {\mu }$-TPIV system.

We image the very midpoint of the array considering a VOI shown in figure 1(b) with the coordinate system triad. Views in three dimensions reference the coloured triad for orientation. For optical access through multiple fluid–solid interfaces throughout the depth of the array, we used RIM between the fluid and solid following a scheme similar to that in Larsson et al. (2018). A continuous 5 mW green laser was passed through the fluid-filled array onto a target 1 m away in a 25 $^\circ$C environment, and the diameter of the distorted beam for different concentrations of glycerol in deionized water was measured. In our configuration, this approach could detect a refractive index (RI) disparity of $10^{-4}$. We obtained a minimum distortion and thus achieved RIM at 89.5 wt% glycerol. A viscoelastic solution was prepared consisting of 130 parts-per-million partially hydrolyzed polyacrylamide (HPAA, $M_W = 18$ MDa, Polysciences Inc., USA), in the RI-matched solvent. For $\mathrm {\mu }$-TPIV measurement purposes, the fluid was seeded with $2\,{\mathrm {\mu }}{\rm m}$ fluorescent particles (excitation/emission 530 nm/607 nm, PS-FluoRed-Particles, Microparticles GmbH, Germany) to a particle concentration of approximately 0.03 particles per pixel at the desired magnification. Ultimately, the RI of the fluid was measured with an Anton-Paar Abbemat MW refractometer operating at 589 nm, yielding RI 1.4582 at 25 $^\circ$C. This is in good agreement with the nominal RI 1.4584 for ${\rm SiO}_2$ at 589 nm (Malitson 1965). Although the refractometer operates at a narrow selection of wavelengths, both materials exhibit normal chromatic dispersion, and the RIM errors incurred from imaging at other wavelengths will partially cancel out (Larsson et al. 2018).

We characterized the shear-rate dependent viscosity and extensional flow characteristics of the HPAA solution using an Anton-Paar MCR stress-controlled rheometer and a Haake capillary breakup extensional rheometer (CaBER), respectively. The fluid is weakly shear-thinning with zero shear viscosity $\eta _0 = 500$ mPa s (figure 2a), and we take a subset of the shear rheology sweep at moderate $\dot {\gamma }$ to fit power-law extrapolations for the first normal stress difference $N_1$ and shear stress $\sigma$ (figure 2b) to evaluate a local $Wi_{L} = N_{1}(\dot {\gamma }_{L})/2\sigma (\dot {\gamma }_{L})$ discussed in § 3. Considering the extensional rheology, the fluid is viscoelastic with relaxation time $\lambda = 1.2$ s (figure 2c), and is strongly strain-hardening with $\eta _{E,app} \approx 340\eta _{0}$ at high strains (figure 2d). Here, the accumulated Hencky strain is $\epsilon _H = 2\ln (d_{0}/d(t))$, with the measured filament diameter $d$, and the initial filament diameter $d_0 = 6$ mm. For flow within the porous medium, the zero shear viscosity $\eta _{0}$ and length scale $w_p/2$ are used to calculate $Re = \rho U_{p}w_{p}/2\eta _{0}$ and the characteristic shear rate $\dot {\gamma } = 2U_{p}/w_{p}$ for $Wi = \lambda \dot {\gamma } = 2\lambda U_{p}/w_{p}$. The Reynolds number $Re \lesssim 10^{-3}$ is considered negligible in the present work. As the array of spheres is cubicly packed, the relevant length scale $L$ for a pitch-based Deborah number is $w_{p}$, thus in the present work, $De = U_{p}\lambda /w_{p} = Wi / 2$.

Figure 2. (a,b) Shear rheology of the HPAA solution, where (b) is at a moderate $\dot {\gamma }$ subset of (a) to fit power-law interpolations for $N_1$ and $\sigma$. (c,d) Extensional rheology of the elastocapillary response for the relaxation time $\lambda$ and strain hardening response.

2.2. Microtomographic PIV ($\mathrm {\mu }$-TPIV)

Tomographic particle image velocimetry (TPIV) is a volumetric flow measurement whereby for a particle-laden flow, numerous camera angles are used to reconstruct a particle volume that is then locally cross-correlated in time to resolve the velocity vector $\boldsymbol {u}$ on a 3-D grid (Elsinga et al. 2006a). It is possible to conduct TPIV at micro-scale ($\mathrm {\mu }$-TPIV) via imaging from the dual LOS afforded by a stereo microscope (see Carlson et al. 2021). Flow was captured by a LaVision FlowMaster $\mathrm {\mu }$-TPIV system in a 25 $^\circ$C environment as time-resolved recordings from dual high-speed cameras (Phantom VEO 410, $1280\times 800$ pixels) at a velocity-dependent frame rate spanning 25–250 Hz such that no particle moved more than 8 pixels between frames when illuminated by a coaxial Nd:YLF laser (527 nm wavelength). Images were pre-processed with an ROI mask spanning 7–5 spheres in $x$–$y$ and filtered with local background subtraction and Gaussian smoothing at $3 \times 3$ pixels. Three-dimensional calibration was achieved by capturing images of a micro-grid in air at the planes $z_{air} = \pm 500$ and $0\,\mathrm {\mu }$m, and mapping the coordinate $z_{air}$ elevation to an equivalent position $z$ in the working fluid due to the RI change. A preliminary coordinate system was achieved by interpolating a third-order polynomial between the calibration planes. The final coordinate system was reached via volumetric self-calibration (Wieneke 2008) of the pre-processed images, which distorts the original grid using the actual particle images to create a refined calibration. By iterating volume self-calibration, we reached a typical maximum disparity of 0.02 voxels. Particle reconstruction was achieved using the Fast MART (multiplicative algebraic reconstruction technique) algorithm implemented in the commercial PIV software DaVis 10.2.1 (LaVision GmbH). This algorithm creates an initial particle volume using multiplicative line of sight (Worth & Nickels 2008; Atkinson & Soria 2009) followed by iterations of Sequential MART (Atkinson & Soria 2009), then sequential motion tracking enhancement (SMTE) for ‘ghost’ suppression (Novara, Batenburg & Scarano 2010; Lynch & Scarano 2015), where spurious ghost particles arise from randomly overlapping LOS and contaminate the cross-correlation (Elsinga, Van Oudheusden & Scarano 2006b). For the present work, we determined that eight iterations of Fast MART resulted in a suitable compromise between vector quality and computation time (6 CPU-minutes per tomogram VOI, $T$), and observed that SMTE converged ghost intensity after a temporal march over 30 volumes. Before temporal cross-correlation, we applied the 3-D masking procedure described below to eliminate the solid region in the reconstruction.

2.3. Tomogram masking

Masking protocols for TPIV range from simple ROI projection to more sophisticated backward projection algorithms such as the visual hull approach (Adhikari & Longmire 2012), but such methods are inapplicable when the RI-matched solid is transparent and the 3-D array has highly overlapping LOS in the back-projection. Therefore we developed an alternative procedure by volume masking directly from the $\mathrm {\mu }$-CT geometry. First, we generate a manifold of the $\mathrm {\mu }$-CT STL about the VOI via a solid Boolean operation in the 3-D modelling software Blender (Blender Online Community 2018), referring to the surface yielded as $S$. Second, we time-averaged a subset of the tomogram volumes $T$ as $\bar {T}$, and further applied $11 \times 11 \times 11$ voxel Gaussian smoothing to generate a reference volume for mask generation. To centre and align $\bar {T}$ with $S$, we leverage the uniform spacing of the array: we project $S$ and $\bar {T}$ onto $x^*$–$z^*$ and $y^*$–$z^*$ planes, and sum about the $x^*$ and $y^*$–$z^*$ axes, respectively. We compare two central $\bar {T}$ maxima against $S$ minima for the $x^*$–$y^*$ axes, and directly align the maxima and minima in $z^*$. An example alignment is shown in figure 3. By weighting multiple minima to locate $S$ in $\bar {T}$ in $x^*$–$y^*$, we obtain sub-voxel mask placement resolution in $x^*$–$y^*$ and single voxel resolution in $z^*$.

Figure 3. Example alignment between the tomogram volume $T$ and the edges of the STL manifold $S$ in the $y^*$–$z^*$ plane. Here, $S$ and $T$ are projected onto the plane, then summed along $y^*$ and $z^*$.

Voxels in the de-warped mask planes of $\bar {T}$ that were found to be inside $S$ via a ray-tracing algorithm are set as masked, and the planes are subsequently warped back into voxel space. We apply the mask planes in DaVis 10.2.1 to $T$ before multi-grid iterative cross-correlation with a final window of $32 \times 32 \times 32$ voxels with 8 voxel steps (75 %) to calculate the velocity field $\boldsymbol {u}$ with components $u$, $v$ and $w$. Finally, we spatiotemporally filtered the vector fields with a second-order polynomial regression across neighbourhoods of $5 \times 5 \times 5 \times 11$ vectors (Scarano & Poelma 2009; Elsinga et al. 2010; Schneiders, Scarano & Elsinga 2017). We can estimate a signal-to-noise ratio (SNR) of this approach by taking the ratio of the average $\bar {T}$ intensity in the fluid (particles and ghosts) and that of the mask region (ghosts and near-boundary particles). In the present work, ${\rm SNR} = 3.0$, indicating that our masked reconstruction is of good quality (${\rm SNR} = 2.2$ without SMTE). Whereas historically, time-resolved PIV has been burdensome technically, modern approaches for TPIV are increasingly reliant on accessing time-resolved quantities such as in SMTE (Beresh 2021).

For a posteriori assessment of our measurement quality, we validated conservation of mass for a time-averaged vector field (Zhang, Tao & Katz 1997). We take flow divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$ assessed relative to an assumption of incompressible flow (relative error $\zeta = (\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u})^2 / \mathrm {tr}(\boldsymbol {\nabla }\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u})$, using the trace $\mathrm {tr}$). At the maximum $Re$, the error is $\zeta = 0.2$ (see figure 4 for joint p.d.f.s of $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$ components and the p.d.f. of $\zeta$), in good agreement with limited LOS TPIV experiments (Kempaiah et al. 2020) and our prior work ($\zeta = 0.25$; Carlson et al. 2021). In addition, a comparable $\zeta$ from our prior work indicates that the fluid–structure RIs are properly matched: the fluid–structure reconstruction results in a vector quality similar to that from a fluid only.

Figure 4. (a) Joint probability density functions (p.d.f.s) of the components of divergence for the solvent at $Re = 10^{-3}$. Points away from the zero divergence plane ($\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} =0$) indicate measurement error, which can be calculated as a relative divergence error parameter $\zeta$. (b) The p.d.f. of $\zeta$.

3. Results and discussion

Pressure and $\mathrm {\mu }$-TPIV measurements were taken over a range of flow rates corresponding to $0.5\leq Wi \leq 8.2$ for the HPAA solution, and $10^{-5}\leq Re \leq 10^{-3}$ for the solvent, considering a VOI midway into the porous array (35 pores). We use the pressure data, which can be interpreted in real time, to gauge the overall flow state and limit resource-intensive $\mathrm {\mu }$-TPIV reconstruction to selected $Wi$ and $Re$. In discussing the flow, we non-dimensionalize the vector field $\boldsymbol {u}$ by the average pore velocity $U_{P}$, and length scales by $w_{p}$ (e.g. $\lvert \boldsymbol {u} \rvert ^*=\lvert \boldsymbol {u} \rvert /U_{p}$, $v^*=v/U_{p}$, $x^* = x/w_{p}$). From the rate-of-strain tensor ${{\boldsymbol{\mathsf{D}}}} = (\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^{\rm T})/2$, we calculate the magnitude $|{{\boldsymbol{\mathsf{D}}}}|$ as $\dot {\gamma }_{L} = \sqrt {2({{\boldsymbol{\mathsf{D}}}}:{{\boldsymbol{\mathsf{D}}}})}$ and reduce it as $\dot {\gamma }^* = \dot {\gamma }_{L} w_p / (2 U_p)$. We introduce the common flow type parameter $\xi$ (Astarita 1979) as a map of the flow topology: $\xi = (|{{\boldsymbol{\mathsf{D}}}}|-|\pmb {\it \varOmega }|)/(|{{\boldsymbol{\mathsf{D}}}}|+|\pmb {\it \varOmega }|)$, where $\pmb {\it \varOmega }= (\boldsymbol {\nabla } \boldsymbol {u}-\boldsymbol {\nabla } \boldsymbol {u}^{\rm T})/2$ is the spin tensor. Thus $\xi = -1$ for flow dominated by solid body rotation, $\xi = 0$ for shear flow, and $\xi = 1$ for extensional flow.

We present the solvent flow field in figure 5 as a baseline case at the maximum $Re$ (${Wi = 0}$). Here, we show $x^*$–$y^*$ slices of $v^*$ and $\dot {\gamma }^*$ as well as numerically advected massless particles (Paraview 5.10; Ahrens, Geveci & Law 2005) coloured by $\lvert \boldsymbol {u} \rvert ^*$ and $\dot {\gamma }^*$ from the background grid. A single ring of particles is injected from $x^* = 0$ at each time step, thus their consistent grouping indicates steady-state flow. The particle fields show fluctuating $\lvert \boldsymbol {u} \rvert ^*$ and $\dot {\gamma }^*$ as they are advected between pore throats and bodies, tracking shear flow through pore throats and extension-dominated flow about the pore bodies consisting of alternating compressional and extensional flow in each pore. Spatially, the fields of $v^*$ and $\dot {\gamma }^*$ are consistent between sequential pores, thus we determine that the VOI is sufficiently deep in the array to consider it fully developed laminar flow.

Figure 5. Instantaneous numerical particle field and $x^*$–$y^*$ slice at $z^* = 0$ for the glycerol solvent flow at $Re = 10^{-3}$. The particle fields are coloured by (a) $\lvert \boldsymbol {u} \rvert ^*$ and (b) the flow type parameter $\xi$, while the slices map the crossflow velocity $v^*$ and rate of strain $\dot {\gamma }^*$, respectively.

3.1. The transition from laminar to elastically turbulent flow

Flow passing through a porous medium at low $Wi$ and low $Re$ is expected to drive a pressure drop that varies approximately linearly with flow rate for the Boger fluids used in our experiments, as per Darcy's law. Thus to delineate our flow regimes, we consider the pressure drop data presented in figure 6. In figure 6(a), we show the time-averaged polymeric contribution to the pressure drop ${\rm \Delta} \bar {p}-{\rm \Delta} \bar {p}_{n}$ between the inlet and outlet of the microchannel. Here, ${\rm \Delta} p$ is the pressure drop (i.e. $p_1 - p_2$) measured for the flow of the polymer solution, and ${\rm \Delta} p_n$ is the pressure drop measured for the Newtonian solvent. Both were recorded over 180 s at each imposed volumetric flow rate $Q$. We see that for low $Wi$ ($0.5\leq Wi \leq 1.8$), the pressure drop obeys the laminar model ($r^2 = 0.99$), with a slight deviation at $Wi = 1.8$, but the polymeric contribution breaks this expectation at higher flow rates $Wi \geq 2.6$, resulting in the often discussed increased flow resistance driven by elastic effects (Marshall & Metzner 1967; Durst et al. 1981; Browne & Datta 2021). We use the departure point of the pressure drop data from the linear trend to estimate the critical $Wi$ as $Wi_{c} = 1.8$, i.e. the highest $Wi$ data point to fall near the trend.

Figure 6. (a) The polymeric contribution to the pressure drop across the flow cell. (b) Power spectra $E_{pp}$ for selected $Wi$ highlighting the gradual strengthening of ET for increasing $Wi$. The dashed line in (a) indicates a linear regression fit below $Wi_{c}$.

To further describe the transition from laminar to unstable flow, we present power spectra $E_{pp}$ of the fluctuating pressure $p' = {\rm \Delta} p - \bar {{\rm \Delta} p}$ against dimensionless frequency $f^* = f\lambda$ in figure 6(b). Despite the pressure drop deviation in figure 6(a) for $1.8\leq Wi \leq 2.6$, when comparing the pressure power spectra of laminar flow at $Wi = 0.5$ to that of an apparently unstable flow at $Wi = 2.6$, one observes that they are quite similar, with slightly increased power in the lower $f^*$ for the unstable state. Evidently, the near-critical unstable flow state is weakly transient (we quantify the flow transience in the Appendix). For increasing $Wi$, the transient state grows in power, and by $Wi = 8.2$, we observe sustained increased power throughout the lower frequency range, with a steep power-law decay exponent $\beta \approx 3$ (which is typical of ET) up to $f^*=1$, i.e. the fluctuation time scale $1/f$ matches that of the polymer relaxation time $\lambda$. We hypothesize that the higher frequency roll-off towards $f^* = 2$ down to the noise floor is due to polydispersity of the HPAA, and generalize that $1/\lambda$ is a common upper bound for capturing the majority of the energy of ET (Browne & Datta 2021). While the power spectra show that a smooth macro-scale transition towards stronger ET occurs with increasing $Wi$, we next consider the local micro-scale evolution.

3.2. Instability: a local curvature balance

In § 1 we introduced a common scaling parameter $M = \sqrt {2\,Wi\,\lambda U/R}$ (McKinley et al. 1996; Pakdel & McKinley 1996) as a generalized predictor for the onset of elastic instability by coupling streamline curvature with normal stresses. The reported literature for $M$ typically uses characteristic values for $Wi$, $U$ and $R$. Instead, we substitute in local parameters such that $M$ can describe the flow on instantaneous streamlines. We first assume that the vector field is quasi-static as sampled (1/250 Hz $\ll 1/f$), then seed streamlines randomly in the $x^*$–$y^*$ plane about $z^* = 0$, and use a Runge–Kutta 4–5 integrator to progress streamlines forwards and backwards in time until they exit the 3-D volume. The streamlines provide a local radius of curvature $R_{L}$ in three dimensions directly, but we linearly interpolate for $U$ and $\dot {\gamma }$ at each streamline from the background grid of $\boldsymbol {u}$ and ${{\boldsymbol{\mathsf{D}}}}$ to provide $\lvert \boldsymbol {u} \rvert$ and the local Weissenberg number $Wi_{L} = N_{1}(\dot {\gamma }_{L})/2\sigma (\dot {\gamma }_{L})$. Thus in the present work, we discuss the local form of the Pakdel–McKinley criterion, $M_{L} = \sqrt {2\,Wi_{L}\,\lambda \,\lvert \boldsymbol {u} \rvert /R_{L}}$.

In figures 7(a–d), we present the outcome of this approach for $M_{L}$ at selected $Wi$ near $Wi_{c}$, in figure 7( f) the p.d.f. of $M_{L}$, and in figure 7(g) $M_{max}$ and $M_{avg}$. For increasing $Wi$ approaching $Wi_{c} = 1.8$, streamline curvature and $M_{L}$ increase, particularly at the throat but also at the midpoint of each pore. Furthermore, the streamlines are pushed outwards incrementally towards the necks between rows for increasing $Wi$. This behaviour is generally symmetric between crossflow pores for $Wi = 0.9, 1.3$. However at $Wi = 1.8$, we see that this symmetry is broken. For any given column of pores in $y^*$, streamlines will tend to flood the neck region asymmetrically from one of the pores, leading to straightened streamlines and reduced $M_{L}$ at the opposing pore.

Figure 7. (a–d) Streamlines of local $M_{L}$ evaluation for selected $Wi$ near $Wi_{c}$, with arrows highlighting asymmetry (sketch in (e)). ( f) The p.d.f. of the streamline fields. (g) The maximum and average $M_{L}$ at each $Wi$. Streamline curvature does increase with $Wi$ above $Wi_{c}$, but in doing so, symmetry is broken as streamlines flood towards neighbouring pores, thereby enforcing stability.

For $Wi = 2.6$ (figure 7d), asymmetric flooding behaviour is shown more distinctly. Regions with this effect have been highlighted with arrows, and a sketch is shown in figure 7(e). Thus we infer that pores with elastically unstable flow will stabilize their peers as the instability manifests as a flooding incursion into neighbouring pores. Notably, this effect also correlates along the streamwise direction, such that rows of pores in $x^*$ show a similar flow response. Above $Wi_{c}$, we can see that local curvature and $M_{L}$ continue to increase with $Wi$, but at a reduced rate (figure 7g), owing to the overall distribution of $M_{L}$ tailing towards higher values but maintaining a similar mean (figure 7g). From $M_{max}$ in figure 7(g), we can estimate our critical $M_{c}$ for this geometry by taking $M_{max}$ at $Wi_c = 1.8$, yielding $M_{c} = 2.5$, with the caveat that although $M_{L}$ is self-consistent for different $Wi$, it will underestimate the theoretical $M$ due to discrete sampling limitations. Nonetheless, $M_{c}$ here is in good agreement with reported values of $M_{c}$ based on $M$ ranging from 1 to 7.

In the transient unstable state, spatial heterogeneity resolves the increased stresses via localized regions of increased flow instability. Although the pressure power spectra indicate that these pockets of instability must be transient, tomogram storage for masking remains computationally burdensome, thus instead we time-resolve the fluctuations of ET at high $Wi$ where stronger low-frequency modes make it feasible to capture multiple state changes in the VOI within a 10 s recording. A sample description of the low-frequency velocity fluctuations is available for $Wi = 1.8$ in the Appendix. Moreover, $Wi_{c} \approx 2$ aligns with $De = \lambda U_{p}/ w_{p} = 1$ by definition, thus transient effects can grow in the array in this regime as polymer chains are advected between pores faster than they can relax.

3.3. Inter-pore heterogeneity

We consider the evolving flow state at $Wi = 8.2$ as a representative case for ET. We support this assertion with selected power spectral densities extracted from a point in the velocity field ($x^* = z^* = 0$, $y^* = -0.75$) to expand upon the prior spectral analysis of the pressure field. In figure 8 we present the global $E_{pp}$ spectra discussed previously, together with local fluctuating crossflow velocity spectra $E_{vp}$, at $Wi = 3.7$ and $Wi = 8.2$. Here, $E_{vp}$ and $E_{pp}$ show the algebraic decays typical of ET: for $E_{pp}(f) \sim f^{-\beta }$ and $E_{vp}(f) \sim f^{-\alpha }$, $\beta \approx 3$ and $\alpha \approx 3.5$, consistent with a recent review (Steinberg 2021). Furthermore, the analysis of Steinberg (2019) deduced a scaling relation $\beta = 2(\alpha -2)$ that supports our experimental data. The Taylor frozen flow hypothesis (Taylor 1938) has been applied reasonably well for parts of ET flow. Specifically, experiments and simulations have noted similarly steep decay exponents for the temporal and spatial power spectra (Burghelea et al. 2005; Burghelea, Segre & Steinberg 2007; van Buel & Stark 2022). Thus our scaling in the frequency domain will also track the wavenumber contribution in the spatial domain, yielding the common observation that ET is governed by the nonlinear interaction of only a few large-scale modes (Burghelea et al. 2005; Steinberg 2019).

Figure 8. Power spectral densities of the fluctuation pressure $E_{pp}(f)$ and the fluctuating crossflow velocity $E_{vp}(f)$ at moderate and high $Wi$. The power law decays for $E_{pp}(f) \sim f^{-\beta }$ and $E_{vp}(f) \sim f^{-\alpha }$ follow a scaling $\beta = 2(\alpha -2)$ (Steinberg 2019).

For a local description of ET flow, in figures 9(a,b) we show streamlines coloured by $\lvert \boldsymbol {u} \rvert ^*$ and $M_{L}$, respectively, to inspect local flow properties between the central rows of pores (labelled as $r_a$ and $r_b$, respectively). Generally, flow through row $r_a$ remains stable from pore columns 0–2, with a flow rate $\lvert \boldsymbol {u} \rvert ^*$ above the volume average, straightened streamlines, and low $M_{L}$. Flow through row $r_b$, however, has an opposite reaction for the same quantities: reduced $\lvert \boldsymbol {u} \rvert ^*$, strongly curving streamlines, and increased local $M_{L}$ across pore columns 1–3, but a reduction in pore 4. Remarkably, we can also see flow through $r_a$ transition to a state similar to flow through $r_b$ by pore 4, in this case driven by out-of-plane streamline curvature towards $+z^*$. Preceding the unstable flow in pore 4, we observe a reduction in $\lvert \boldsymbol {u} \rvert ^*$ in pore 3 as either the elastic instability or other rheological parameters drive increased pore resistance. Finally, we also see the same correlation behaviours noted in the prior section: streamwise correlation and crossflow inverse correlation of pores with unstable flow. The inversely correlated state is particularly compelling, with no $y^*$-adjacent pair of pores responding with mirrored unstable flow. This further supports our hypothesis that pore flooding can promote stability in adjacent pores. In pores farther from the unstable region, however, we do observe adjacent stable flow states, such as at $x^* = y^* = 0$ in figure 9. Therefore, the adjacent peer stabilization effect may be alternatively classed as a crossflow exclusion of instability. Streamwise correlation generally holds true, apart from the transition states noted in pores 3 and 4 as $r_a$ goes unstable and $r_b$ reduces curvature, respectively. A likely explanation for the streamwise consistency is flow history, which we will investigate by continuing this analysis throughout time in § 3.4.

Figure 9. Instantaneous streamlines for the HPAA flow under a state of ET at $Wi = 8.2$, coloured by (a) the dimensionless velocity magnitude $\lvert \boldsymbol {u} \rvert ^*$, and (b) the Pakdel–McKinley criterion $M_{L}$. We have numbered columns of pores along $x^*$.

To describe the laminar and unstable flow states, we take a snapshot also at $Wi = 8.2$, but at an instant when $r_a$ and $r_b$ are in consistent states along $x^*$. Figure 10 presents both the local $M_{L}$ streamlines previously described together with planes of $\lvert \boldsymbol {u} \rvert ^*$, $v^*$, $\xi$ and $\dot {\gamma }$. Here, $r_a$ and $r_b$ have almost fully reversed their states from the figure 9 snapshot, with $r_a$ uniformly unstable between pores along $x^*$ and in $r_b$ repeating pores with laminar flow. Hence we can relate the instability to local rheological effects. As discussed previously, pores with unstable flow lead consistently to increased flow resistance (locally decelerated) compared to the laminar channel, but importantly, we see here that the laminar channel $r_b$ is also shear-dominated with strong rate-of-strain. The shear-thinning behaviour of the fluid is likely in effect, leading to a locally accelerated flow reaching $\lvert \boldsymbol {u} \rvert ^* = 3.6$. Along the unstable zone, we observe competing influences: while the fluid is strongly extensionally strain hardening, and these pores manifest with extensionally dominant flow, the rate-of-strain is locally reduced. Thus while strain hardening may be an influence in the evolution of the local instability, we cannot confirm this hypothesis in the present work.

Figure 10. (a,b) Instantaneous streamlines for the HPAA flow under a state of ET at $Wi = 8.2$, together with (c–f) $x^*$–$y^*$ slices showing $\lvert \boldsymbol {u} \rvert ^*$, $\xi$, $v^*$ and $\dot {\gamma }$, respectively.

3.4. Feedback between streamline curvature and flow resistance

To describe the temporal evolution of the unstable flow state, we first refer back to the fluctuating pressure power spectra from figure 6(b). As we mentioned previously, $f^* \approx 1$ appears to be a natural upper limit for the frequency of the spatiotemporal modes of ET as fluctuations approach $1/\lambda$ but decay rapidly at higher frequencies (figure 8). Therefore, we compare snapshots with a temporal stride of approximately $\lambda$, but with some variation permitted to capture disturbances evolving away from $\lambda$.

Selected streamline fields are presented in figure 11 to describe generally instability propagation forward in time. We start the sequence in figure 11(a), where an individual pore in the VOI is unstable (3, the last pore). All of the pores upstream of this point are stable in this row. Forward in time at $t^* = 0.7$ (figure 11b), we see that pores 2 and 3 are now unstable; the unstable region has evidently propagated upstream. This trend continues for $t^* = 1.8$ and $2.8$ (figures 11c,d), with the full row now unstable at $t^* = 2.8$. Upstream propagation continues at $t^* = 3.9$ (figure 11e), with relaminarization occurring in pore 3. Clearly, flow history is important in this context, for both the onset of instability as well as the termination thereof. We support this notion by considering the streamline curvature of pores upstream of the local maximum $M_{L}$. The transition from a laminar state to a strongly unstable flow is spatially smooth, with streamline curvature and elastic stresses building for pores progressing along increasing $x^*$ towards a pore with a locally maximum state of instability (protruding streamlines). Positive feedback between streamline curvature and flow resistance is a likely explanation: in a given pore, local streamlines respond with increased curvature due to raised resistance in neighbouring pores downstream. In a given pore with unstable flow, streamlines protrude from the pore centre (figure 7e), and this excursion becomes more pronounced over time while being accompanied by locally increased flow resistance. Thereafter, flow in an adjacent pore upstream responds to the raised resistance with streamline protrusion, hence the instability propagates upstream in the array.

Figure 11. Snapshot streamline fields of $M_{L}$ over selected $t^*$ at $Wi = 8.2$. Note that the black line in (a) indicates $y^* = -0.75$, which is used in figure 12.

To further explain the instability propagation, we present space–time diagrams in figure 12 extracted along $y^* = \pm 0.75$ from separate recordings at the same $Wi = 8.2$. For reference, we also show the advection length scale $U_{p}\lambda /w_{p}$ overlaid on the figure. Here, we see that for a given point $x^*$, flow maintains a consistent laminar or unstable state over long $t^*$. Whereas for a time $t^*$, the fluctuating state persists over $x^*$ with a scale on a par with the advection length scale. Polymer chains in the flow would thus be in a state where strain is accumulated throughout a cohesive pocket of instability. Crossflow fluctuation may propagate rapidly upstream, as it does near $t^{*}_{0} = 2$ for $x^* = 2\unicode{x2013}3$, or the laminar state could linger as it does along $t^{*}_{0} = 2\unicode{x2013}4$ for $x^* = 1$. Similar behaviour is noted in the termination of instability shown in figure 12(b), where, for example, fluctuations decay along low $x^* = 0\unicode{x2013}2$ over a long time scale surpassing the recording limit. These time scales are not surprising when we compare the space–time diagrams to the fluctuating pressure power spectra from figure 6, where the lowest frequencies continued to gain power for decreasing $f^*$ towards the lower limit $0.06 f\lambda$ (period $17\lambda$). Thus although an upstream propagation speed of instability can be estimated roughly as $\approx$$-0.3U_p$, the progression from laminar to unstable flow is spatially smooth but not temporally, in support of the spectral characteristics of ET.

Figure 12. Space–time diagrams along $y^* = \pm 0.75$ of the crossflow velocity magnitude ratio $|v|/|\boldsymbol {u}|$ from separate recordings at $Wi = 8.2$. We also include the advection length scale $U_{p}\lambda /w_{p}$ as the pink bar.

From our simplified ordered porous medium, we have elucidated viscoelastic flow effects relevant to porous media flow in particular, and to the evolution of vessel-scale modes for ET in general. Namely, we have highlighted the importance of spatial correlation via flow history effects as well as via communication between pore flow states. An area for further research, however, would be probing longer time scales; such are not feasible for a volumetric time-resolved study given current processing times, but a 2-D study of ET in our ordered porous medium would be a promising direction for future work, to better quantify the propagation nature of unstable flow near $Wi_{c}$. Furthermore, while ET has been quantified in 3-D randomly packed media (Browne & Datta 2021), packing structure and disorder are known to be significant for the onset of instability for viscoelastic flow (Walkama et al. 2020; Haward et al. 2021) and for chaotic mixing for creeping porous media flow (Turuban et al. 2018, 2019), and these are thus topics of ongoing research.

4. Conclusions

From time-resolved measurements of dilute polymer flow through a refractive- index-matched novel arrayed porous medium, we have produced the first spatiotemporally complete description of elastic turbulence (ET). By reconstructing flow over multiple pores simultaneously, we can capture the local evolution of viscoelastic instability as flow evolves over a time span surpassing that of the polymer relaxation time $\lambda$. Above a critical $Wi_{c} = 1.8$, a purely elastic instability occurs, leading to increased resistance, strongly curving streamlines and localized pore flooding, the last of which is an effect that actually stabilizes flow in neighbouring pores in crossflow by causing a reduction in streamline curvature there. The strongly fluctuating flow is accompanied by increased drag and in the spectral domain, steep power-law decays for velocity and pressure matching proposed ET scaling for $\alpha$ and $\beta$. In this state, we describe unsteady yet coexisting regions of laminar flow and unstable flow. We relate the spatially heterogeneous laminar zones to the local shear-thinning behaviour of the solution, but see low rate-of-strain with extensional flow in the resisting pores. We hypothesize that the pockets of instability are transient in nature and driven by locally raised downstream resistance acting as positive feedback on upstream streamline curvature. Whereas ET is typically characterized at the vessel scale, these new findings considering a flow volume developing over several vessel scales elucidate the significance of spatiotemporal correlation on the local evolution of ET.