1. Introduction
The flow of polymer solutions through porous media is widely used across many applications, ranging from groundwater remediation and enhanced oil recovery to filtration and polymer processing (Boudreaux & Cuculo Reference Boudreaux and Cuculo1976; Bourgeat, Gipouloux & Marusic-Paloka Reference Bourgeat, Gipouloux and Marusic-Paloka2003; Smith et al. Reference Smith, Silva, Munakata-Marr and McCray2008; Sorbie Reference Sorbie2013; Pogaku et al. Reference Pogaku, Mohd Fuat, Sakar, Cha, Musa, Awang Tajudin and Morris2018; Firozjaii & Saghafi Reference Firozjaii and Saghafi2020). In these applications, a polymeric fluid must navigate tortuous flow paths set by the geometry of the pore space. If the flow is sufficiently fast, the constituent polymer molecules become stretched as they are transported along the curved streamlines, and the elastic stresses thereby generated develop and give rise to a flow instability (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990; Larson Reference Larson1992; Shaqfeh Reference Shaqfeh1996; McKinley, Pakdel & Öztekin Reference McKinley, Pakdel and Öztekin1996; Pakdel & McKinley Reference Pakdel and McKinley1996; Pathak, Ross & Migler Reference Pathak, Ross and Migler2004; Browne et al. Reference Browne, Shih and Datta2020b
; Steinberg Reference Steinberg2021; Datta et al. Reference Datta2022; Kumar et al. Reference Kumar, Guasto and Ardekani2022a
). The Reynolds number
$\textit{Re}$
associated with such flows is typically
$\ll 1$
(Browne et al. Reference Browne, Shih and Datta2020b
), and so, this instability is known as a purely elastic instability. Its onset can instead be parametrised using the Weissenberg number,
$\textit{Wi} = N_1/ (2\sigma )$
, which compares the strength of elastic stresses, described by the first normal stress difference
$N_1$
, to viscous stresses, described by the shear stress
$\sigma$
.
What are the characteristics of elastic instabilities in a porous medium? Despite its apparent simplicity, the answer to this question is unclear. One reason is because the flow is highly sensitive to pore space geometry: studies in idealised one- and two-dimensional (1-D and 2-D) microfluidics (Anbari et al. Reference Anbari, Chien, Datta, Deng, Weitz and Fan2018; Browne et al. Reference Browne, Shih and Datta2020b ; Datta et al. Reference Datta2022) using a variety of polymeric fluids have uncovered a dizzying array of different flow behaviours in different geometries, including symmetry-breaking (Arratia et al. Reference Arratia, Thomas, Diorio and Gollub2006; Dubash, Cheung & Shen Reference Dubash, Cheung and Shen2012; Haward & McKinley Reference Haward and McKinley2013; Haward, Hopkins & Shen Reference Haward, Hopkins and Shen2020; Qin et al. Reference Qin, Ran, Salipante, Hudson and Arratia2020; Hopkins et al. Reference Hopkins, Shen and Haward2022b ), quasi-periodicity (Zhao, Shen & Haward Reference Zhao, Shen and Haward2016; Hopkins et al. Reference Hopkins, Haward and Shen2022a ), the formation of pronounced eddies and wakes (Shi & Christopher Reference Shi and Christopher2016; Zhao et al. Reference Zhao, Shen and Haward2016; Kawale et al. Reference Kawale, Marques, Zitha, Kreutzer, Rossen and Boukany2017; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019b ; Haward et al. Reference Haward, Hopkins and Shen2020; Walkama, Waisbord & Guasto Reference Walkama, Waisbord and Guasto2020; Browne et al. Reference Browne, Shih and Datta2020a ; Kumar et al. Reference Kumar, Aramideh, Browne, Datta and Ardekani2021; Haward, Hopkins & Shen Reference Haward, Hopkins and Shen2021; Hopkins et al. Reference Hopkins, Haward and Shen2022a , Reference Hopkins, Shen and Hawardb ; Dzanic et al. Reference Dzanic, From, Gupta, Xie and Sauret2023; Chen & Datta Reference Chen and Datta2025), and chaotic dynamics reminiscent of turbulence (Groisman & Steinberg Reference Groisman and Steinberg2000, Reference Pathak, Ross and Migler2004; Zilz et al. Reference Zilz, Poole, Alves, Bartolo, Levaché and Lindner2012; Casanellas et al. Reference Casanellas, Alves, Poole, Lerouge and Lindner2016; Qin & Arratia Reference Qin and Arratia2017; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019a ; van Buel & Stark Reference van Buel and Stark2022). This broad diversity of ways in which the flow instability manifests across different geometries can be challenging to make sense of. As a result, studies probing the unsteady flow in porous media with defined and systematically tunable geometries are critically needed. Addressing this gap in knowledge is not just fundamentally interesting, but also practically important: there is tremendous interest in harnessing this elastic instability to homogenise fluid flow (Browne et al. Reference Browne, Huang, Zheng and Datta2023), enhance solute mixing and dispersion (Groisman & Steinberg Reference Groisman and Steinberg2001; Burghelea et al. Reference Pathak, Ross and Migler2004; Pathak et al. Reference Pathak, Ross and Migler2004; Scholz et al. Reference Scholz, Wirner, Gomez-Solano and Bechinger2014; Jacob, Krishnamoorti & Conrad Reference Jacob, Krishnamoorti and Conrad2017; Ducloué et al. Reference Ducloué, Casanellas, Haward, Poole, Alves, Lerouge, Shen and Lindner2019; Browne & Datta Reference Browne and Datta2024), and mobilise trapped immiscible fluids (Clarke et al. Reference Clarke, Howe, Mitchell, Staniland, Hawkes and Leeper2015, Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016; De et al. Reference De, Krishnan, van der Schaaf, Kuipers, Peters and Padding2018; Shakeri, Jung & Seemann Reference Shakeri, Jung and Seemann2021) in porous media. Furthermore, the onset of instability at the pore scale has been linked to a substantial increase in flow resistance at large scales (Arora, Sureshkumar & Khomami Reference Arora, Sureshkumar and Khomami2002; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland, Hawkes and Leeper2015; Browne & Datta Reference Browne and Datta2021) and, thus, has major implications for macroscopic operating conditions. However, in all these cases, a straightforward link between imposed flow conditions, pore space geometry and the characteristics of the flow instability remains lacking.
A key conceptual advance arising from prior 2-D studies is the importance of stagnation points – points at which the local fluid velocity is zero that arise as streamlines diverge/converge around solid boundaries. Stagnation points are therefore dominated by locally extensile flow that generates large tensile stresses in the polymeric fluid that control the onset and strength of the elastic instability (Öztekin et al. Reference Öztekin, Alakus and McKinley1997; Arratia et al. Reference Arratia, Thomas, Diorio and Gollub2006; Poole, Alves & Oliveira Reference Poole, Alves and Oliveira2007; Shi & Christopher Reference Shi and Christopher2016; Zhao et al. Reference Zhao, Shen and Haward2016; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019b ; Haward et al. Reference Haward, Kitajima, Toda-Peters, Takahashi and Shen2019, Reference Haward, Hopkins and Shen2020, Reference Haward, Hopkins and Shen2021; Hopkins et al. Reference Hopkins, Haward and Shen2022a , Reference Hopkins, Shen and Hawardb ). The density and location of these points is dictated by boundary geometry. Hence, identifying when and where stagnation points arise in media of different geometries may provide a straightforward way to disentangle the influence of pore space geometry on the unsteady flow of a polymer solution. Though studies in 2-D have provided useful insights, pore space topology – and thus, the occurrence of stagnation points – is fundamentally different in a three-dimensional (3-D) porous medium (Lester, Dentz & Le Borgne Reference Lester, Dentz and Le Borgne2016; Marafini et al. Reference Marafini, La Rocca, Fiori, Battiato and Prestininzi2020). Therefore, we ask: Where do stagnation points arise in 3-D porous media of different geometries and what are the implications for the unsteady flow of a polymeric fluid?
Stagnation points in flow past spheres. (a) Flow past an isolated sphere with polar stagnation points marked by red
$\times$
s. (b) Flow past consolidated spheres joined together by a neck, producing new stagnation points marked by yellow
$\times$
s. (c) Image of a consolidated simple cubic packing, taken from top-down in air.
Here, we address this question by directly visualising the pore-scale flow of an elastic polymer solution inside 3-D porous media of defined geometries. We use microfabricated simple cubic (SC) and body-centred cuboid (BC) packings of monodisperse glass spherical grains, which circumvent the geometric complexity of disordered media in which the features of the flow instability can differ from pore to pore (Browne & Datta Reference Browne and Datta2021). Stokes flow around a single spherical grain has two stagnation points: one at the upstream pole, where the fluid streamlines diverge around the grain, and the other immediately downstream, where the streamlines reconverge (figure 1 a). These points would then generate strong local extensional flow that can stretch polymers and thereby control the elastic instability (Haward & Odell Reference Haward and Odell2004). However, in certain ordered 3-D packings, which are often conceptualised as having point-like contacts between adjacent grains, these polar stagnation points are screened from the flow due to the presence of the contacting upstream and downstream grains. Rather, in practice, the contacts between adjacent grains are not discrete points, but have extended catenoid-like geometries arising from grains fusing together (figures 1 b and 1 c) – as is the case during consolidation and diagenesis, which are processes inherent in most natural and engineered media. We find that stagnation points are generated at these contacts instead, and these new stagnation points are what control the onset and features of the elastic flow instability, which manifest as a transition from steady to unsteady, time-dependent flow. Moreover, by combining our flow visualisation with pressure drop measurements, we show that the extensional flow associated with these contact-associated stagnation points causes the macroscopic flow resistance to sharply increase, helping to rationalise previous pressure drop measurements (Haward & Odell Reference Haward and Odell2003; Odell & Haward Reference Odell and Haward2006). Our work thus reveals the previously overlooked, pivotal role of grain contacts in shaping elastic flow instabilities in 3-D porous media.
2. Materials and methods
We investigate flow through two different ordered packings of spherical grains (diameter
$D_{\!p} = 800 \: \unicode {x03BC} \mathrm{m}$
) with an overall square cross-section of area
$A_c=16 \: \mathrm{mm^2}$
, as schematised in figure 2(a). The packings are made by selective-laser induced etching of an open-faced slab of fused silica, as detailed previously (Carlson et al. Reference Carlson, Toda-Peters, Shen and Haward2022). They have an SC or BC unit cell; the latter has a rectangular cuboid unit cell, with its longest edge along the direction of imposed flow, instead of the uniform cubic cell typical of an ideal BC packing, to maintain the same number of unit cells in each cross-section. The SC array contains 5
$\times$
5
$\times$
47 packed grains, giving a total medium length of
$\Delta L = 3.84$
cm and porosity
$\phi = 0.48$
. The BC array contains 5
$\times$
5
$\times$
23 grains (excluding the sphere in the centre of each unit cell), giving a total medium length of
$\Delta L= 2.74$
cm and porosity
$\phi =0.28$
. Importantly, the contacts between adjacent grains are not discrete points; instead, due to the microfabrication process, they are small ‘necks’ of diameter
$D_{\mathit{neck}} = 170 \: \unicode {x03BC} \mathrm{m}$
, measured from confocal imaging, which are representative of granular contacts in many natural and engineered media. To enable polymeric fluid flow through the pore space, we glue flexible Tygon tubing and attach plastic connectors (McMaster-Carr) to each microfluidic device.
Experimental set-up and fluid rheology. (a) Schematic of experimental set-up. The microfluidic device containing the microfabricated packing (not to scale) is mounted on the stage of a confocal microscope. A pressure transducer measures the bulk pressure drop across the entire medium. (b) Steady-shear rheology of the 300 ppm HPAM solution. Linear fits correspond to power law models for both the shear stress,
$\sigma$
, and first normal stress difference,
$N_1$
. Error bars represent one standard deviation from averaging over three replicate flow curves, using new fluid samples for each replicate. In some cases, error bars may be smaller than symbol size. The range of
$N_1$
data is limited due to measurement noise of the normal force at lower shear rates. (c) Sample stress relaxation measurement. A controlled shear rate of
$1 \: \mathrm{s}^{-1}$
is applied for 60 s and immediately stopped. The longest polymer relaxation time
$\lambda$
is found by fitting an exponential stress relaxation, assuming a single relaxation mode.
While such 3-D media are typically opaque, we build on our prior work and formulate a polymeric fluid whose refractive index is matched to that of the silica; thus, each medium becomes transparent when saturated with the fluid, enabling us to visualise the motion of fluorescent tracers inside the 3-D pore space and thereby map the flow field in situ using particle image velocimetry (PIV) (Datta et al. Reference Datta, Chiang, Ramakrishnan and Weitz2013; Browne & Datta Reference Browne and Datta2021; Carlson et al. Reference Carlson, Toda-Peters, Shen and Haward2022). We prepare the polymeric fluid by first dissolving partially hydrolysed polyacrylamide (HPAM; 18 MDa; 30 % carboxylated monomers; Polysciences) in ultrapure Millipore water, and adding glycerol (Sigma–Aldrich), dimethyl sulfoxide (DMSO; Sigma–Aldrich) and sodium chloride (NaCl; Sigma–Aldrich). The final solution composition contains 300 ppm HPAM, 82.6 wt % glycerol, 10.4 wt % DMSO and 1 wt % NaCl. This solution falls in the dilute concentration regime,
$c/c^* = 0.5$
, where
$c^* = 600 \pm 300$
ppm is the polymer overlap concentration (Browne & Datta Reference Browne and Datta2021).
To prepare each porous medium before performing flow visualisation, we sequentially flush the device with milliQ water, isopropanol (IPA; Fisher Chemical) and milliQ water again, before filling the medium with glycerol. We then displace the glycerol with the polymer solution at a constant flow rate
$Q = 3 \: \mathrm{mL\,hr}^{-1}$
using a Harvard Apparatus PHD 2000 syringe pump for at least 3 hours to pre-saturate the medium with the polymer solution. During each flow visualisation experiment, we apply a ramp of increasing flow rate, where we apply flow for at least 90 minutes at each new flow rate before beginning imaging to achieve an equilibrated flow state. Due to mechanical pumping limitations, we are limited to examining flow conditions corresponding to
$1\lt \textit{Wi}\lt 5$
in our experiments. To visualise the flow, we seed the fluid with 7 ppm of fluorescent 1
$\unicode {x03BC}$
m carboxylate-modified polystyrene tracer particles (Invitrogen), which have excitation at 480–510 nm and emission at 515 nm. The carboxylate functionalisation helps mitigate deposition onto the silica grain surfaces and the particle concentration was additionally chosen to minimise deposition while maintaining PIV accuracy using the selected optical imaging settings. Particles are excited by a 488 nm laser and detected by a 500–550 nm sensor. For imaging, the medium is placed on the stage of a Nikon A1R+ laser scanning confocal fluorescence microscope. We use a 10
$\times$
objective to interrogate a 1272.79
$\times$
1272.79
$\unicode {x03BC}$
m field of view at 512
$\times$
512 pixels (2.46
$\unicode {x03BC}$
m per pixel). The optical sectioning provided by the confocal microscope and fluid refractive index matching enables clear imaging of a
$z$
-slice of the flow field within the 3-D medium. We use the confocal resonant scanner to capture videos at 30 fps with an optical section thickness of 7.39
$\unicode {x03BC}$
m at depths
$400 \lt z \lt 1200 \: \unicode {x03BC}$
m from the bottom of the medium. We image 10 randomly selected pores across the porous medium and do not observe significant differences in flow behaviour for pores beyond
${\sim} 1\ \mathrm{to}\ 2$
sphere layers near the inlet, presumably due to the identical geometry of each pore throughout the packing. Thus, we select the pore located in the centre of the medium for the majority of the subsequent analysis presented in the results.
We use a built-in MATLAB particle image velocimetry toolbox (PIVlab) to obtain 2-D temporally resolved velocity fields (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014). We use the fast Fourier transform (FFT) window deformation PIV algorithm with four passes of decreasing interrogation area (64–32–16–16 pixels) and 50 % overlap to obtain velocity fields,
$\boldsymbol{u}=(u,v)$
, on a spatial grid,
$\boldsymbol{x}=(x,y)$
, with spatial resolution of
$\Delta x = \Delta y = 20.2 \: \unicode {x03BC}$
m and a temporal resolution of 33 ms. We use image processing software (ImageJ) to generate fluid pathline videos and time-averaged pathline images by applying a moving average of mean pixel intensity over a selected number of frames.
To characterise the macroscopic flow resistance, we measure the pressure drop across the entire medium by connecting a differential pressure transducer (0–30 psi, Omega PX26) in parallel to the device. We apply a constant flow rate using a syringe pump and wait 60 minutes before recording the pressure signal while sampling at a rate of 10 Hz for 90 minutes at each applied flow rate. We calculate a time-averaged pressure signal, adjusted for measurement noise by subtracting off the average pressure signal measured under no-flow conditions. We then determine an apparent viscosity from the pressure drop measurements using Darcy’s law:
$\eta _{\mathit{app}}=\Delta Pk/ ((Q/A_c)\Delta L )$
, where
$k$
is the absolute permeability of the porous medium. We determine the permeability using a flow rate ramp of pure water in each medium, extracting
$k$
from the slope of the linear fit to the pressure,
$\Delta P$
, versus flow rate,
$Q$
. The permeabilities of each medium are
$k_{\mathit{SC}}=567 \pm 9 \: \unicode {x03BC}\text{m}^2$
and
$k_{\mathit{BC}}=197 \pm 2 \: \unicode {x03BC} \mathrm{m}^2$
.
We characterise the bulk solution rheology using an Anton-Paar MCR501 rheometer fitted with a truncated cone-plate measuring tool (CP50-2: 50 mm, 2
$^\circ$
, 53
$\unicode {x03BC} \mathrm{m}$
gap) operating at a controlled temperature of 25
$^\circ$
C. Figure 2(b) presents steady-shear rheology measurements of the shear stress (
$\sigma$
) and first normal stress difference (
$N_1$
) over shear rates
$0.1 \lt \dot {\gamma } \lt 50 \ \mathrm{s^{-1}}$
. As shown by the circles, the shear stress increases approximately linearly with shear rate; a power law fit (
$\sigma \sim \dot {\gamma }^{n_{\sigma }}$
) has an exponent
$n_{\sigma } = 0.93$
, indicating that the polymer solution is very weakly shear-thinning, so we do not expect effects of shear-thinning or shear-banding to influence the flow behaviour. We use the shear stress at the lowest measured shear rate (
$\dot {\gamma }=0.1 \: \mathrm{s^{-1}}$
) to estimate the zero-shear viscosity as
$\eta _0 \approx \sigma /\dot {\gamma }|_{0.1 \: \mathrm{s^{-1}}}= 0.4 \: \mathrm{Pa \: s}$
. The Newtonian solvent viscosity is
$\eta _s=0.13 \: \mathrm{Pa \: s}$
. The first normal stress difference also increases with shear rate as a power law (
$N_1 \sim \dot {\gamma }^{n_{N_1}}$
) with exponent
$n_{N_1} = 1.38$
, indicating that the polymer solution is elastic. We also use shear stress relaxation measurements to estimate the solution relaxation time
$\lambda$
(Liu et al. Reference Liu, Jun and Steinberg2007, Reference Liu, Jun and Steinberg2009); an example is shown in figure 2(c). Averaging over three replicates and fitting the data to a single-mode exponential decay (
$\sigma \sim e^{-t/\lambda }$
) yields a longest relaxation time of
$\lambda = 0.4 \pm 0.1$
s.
Throughout this paper, we use the term ‘elastic instability’ to describe a purely elastic flow instability occurring at negligible Reynolds number. Here, we focus on identifying how spatially extended grain contacts influence the onset of unsteady flow, irrespective of the precise nature of the instability. Identifying the nature of this instability and exactly how it produces unsteady flow will be a useful direction for future research.
3. Results
3.1. Pore-scale imaging reveals the onset of an elastic instability above a threshold
$\textit{Wi}$
First, we use in situ visualisation to probe how the pore-scale polymer solution flow changes with increasing flow rate in both ordered packings. The characteristic interstitial shear rate is given by
$\dot \gamma _I=U_I/\sqrt {k/\phi }$
, where
$U_I=Q/ (\phi A_c )$
and
$\sqrt {k/\phi }$
are the characteristic interstitial flow speed and pore length scale, respectively (Zami-Pierre et al. Reference Zami-Pierre, de Loubens, Quintard and Davit2016; Berg & van Wunnik Reference Berg and van Wunnik2017). We estimate the Reynolds number as
$\textit{Re}=\rho U_I D_{\!p}/\eta _I$
, where
$\rho$
is the fluid density and
$\eta _I$
is the bulk shear viscosity evaluated at the interstitial shear rate
$\dot \gamma _I$
. Here, we use the grain diameter
$D_{\!p}$
as the characteristic length scale to provide an upper bound estimate of
$\textit{Re}$
. We also compute the Weissenberg number,
$\textit{Wi}=N_1/(2\sigma )$
, directly from the bulk shear measurements of
$N_1$
and
$\sigma$
evaluated at
$\dot \gamma _I$
(figure 2
b). Our experiments probe
$10^{-5} \lt \textit{Re} \lt 10^{-3} \ll 1$
and
$1 \lt \textit{Wi} \lt 5$
.
Elastic instability in SC packing: 3-D fluctuating eddies. (a) Schematics contrasting an idealised SC unit cell (top) containing pure point contacts with the consolidated SC unit cell (bottom) used in this work. The small region outlined in red between spheres marks an example of a neck-like inter-grain contact. Yellow
$\times$
s demarcate stagnation points arising at the contacts. Pink
$\times$
s also mark contact-generated stagnation points, but for contacts that are rotated 90
$^\circ$
with respect to the
$z$
-plane shown in the pathline images. (b) Pathline images corresponding to the
$z$
-plane containing the inter-grain contacts (example contact circled in red). (i) Steady flow field at
$\textit{Wi} = 1.5$
. Pathlines are generated by averaging successive frames over 20 s intervals. Scale bar is 200
$\unicode {x03BC}$
m. (ii) Unsteady flow field at
$\textit{Wi} = 4.5$
for three different time points. Highly fluctuating eddies develop between adjacent contacts. Pathlines are generated by averaging successive frames over 2 s intervals. White arrows indicate cross-flow ‘burst’ events between adjacent pores. (c) Pathline images in the
$z$
-plane located 1/2-sphere diameter above the
$z$
-plane in panel (b). (i) Steady flow field at
$\textit{Wi} = 1.5$
. (ii) Unsteady flow field at
$\textit{Wi} = 4.5$
for three different time points, where the flow field fluctuates chaotically.
3.1.1. Simple cubic packing: 3-D fluctuating eddies
The idealised representation of the SC unit cell assumes point contacts between adjacent grains, as schematised in the top of figure 3(a). By contrast, real-world media typically have spatially extended contacts between grains due to consolidation or diagenesis, in which grains fuse together due to the weight of overlying material or biogeochemical processes, respectively; indeed, a 3-D packing with purely point contacts would not sustain rigidity under flow (Deresiewicz Reference Deresiewicz1958; Biot Reference Biot1963; Winkler Reference Winkler1983; Surdam et al. Reference Surdam, Crossey, Hagen and Heasler1989; Mueth, Jaeger & Nagel Reference Mueth, Jaeger and Nagel1998; Makse et al. Reference Makse, Gland, Johnson and Schwartz2004; Burley & Worden Reference Burley and Worden2009; Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013; Horabik & Molenda Reference Horabik and Molenda2016). Our microfabrication process recreates this inherent feature of granular packings. In particular, the unit cell of our microfabricated SC packing has neck-like contacts between adjacent grains – as shown by the red outlines in the schematic at the bottom of figure 3(a) and confocal micrographs in figure 3(b). As we detail in the following, this simple geometric feature has pivotal implications for the unsteady flow behaviour.
At low
$\textit{Wi} = 1.5$
, the flow in the interstitial pore space is steady, as shown by the steady fluid pathlines shown in figure 3(b-i). By contrast, for large enough
$\textit{Wi}$
, an elastic instability is generated by the additional stagnation points directly upstream and downstream of the inter-grain contacts (yellow
$\times$
in figures 3
a and 3
b). Eddies form on the upstream faces of the contacts and grow in length with increasing
$\textit{Wi}$
, localising in the regions between neighbouring contacts, as shown in supplementary movie 1 available at https://doi.org/10.1017/jfm.2026.11458 for a range of
$1.5\lt \textit{Wi}\lt 4.5$
. Above a threshold
$\textit{Wi}=\textit{Wi}_c\approx 3$
, these eddies become highly unsteady, as exemplified by the three different panels in figure 3(b-ii) for the case of
$\textit{Wi} = 4.5$
. They exhibit large fluctuations in size due to ‘bursts’ of transverse cross-flow between pores, as indicated by the white arrows, corroborating the previous observations of Carlson et al. (Reference Carlson, Toda-Peters, Shen and Haward2022).
The variation of this flow field in 3-D can be seen by imaging other
$z$
planes of the SC unit cell. For example, figure 3(c) shows additional pathline images of the flow for a
$z$
-plane located half a grain diameter above that shown in figure 3(b); imaging this plane is equivalent to imaging through an eddy due to the symmetry of the packing. Stagnation points generated at inter-grain contacts rotated 90
$^\circ$
relative to the plane of imaging are denoted by the pink
$\times$
. Again, at sufficiently large
$\textit{Wi}$
, the eddies are highly unsteady, as shown by the crossing pathlines and distinct changes in mean flow direction at different time points in figure 3(c-ii). Supplementary movies 2 and 3 show the flow behaviour at
$\textit{Wi}=1.5$
(steady) and
$\textit{Wi} = 4.5$
(unsteady) for additional
$z$
-planes spaced approximately 50
$\unicode {x03BC}$
m apart that span half of the characteristic unit cell volume. The pathline images and supplementary movies show evidence of tracer particle deposition on grain surfaces; however, flow resistance measurements of tracer and tracer-free polymer solution in both packing geometries show identical resistance behaviour and onset
$\textit{Wi}_c$
. Thus, we do not expect any surface roughness effects resulting from tracer deposition to influence the onset of instability in our experiments.
Spatiotemporal features of the SC eddy instability. (a) Left: Normalised r.m.s. velocity fluctuation maps for steady (
$\textit{Wi}=1.5$
, top left) and unsteady (
$\textit{Wi}=4.5$
, bottom left) flow states. Scale bar is
$200 \: \unicode {x03BC}$
m. Red circles mark examples of inter-grain contacts. Right: Kymographs of the normalised velocity fluctuations corresponding to the column of pixels boxed in the colourmaps on the left. Time
$t$
is reported as the number of injected pore volumes,
$\tau _{PV}=t/(\mathit{PV}/Q)$
, where one pore volume (PV) is defined as the void volume of one unit cell of the packing. (b) Normalised r.m.s. velocity fluctuation maps for additional
$z$
-planes at
$\textit{Wi}=4.5$
spanning half of the unit cell. Sequential
$z$
-planes are spaced approximately
$50 \: \unicode {x03BC}$
m apart.
To more clearly represent these complex dynamics, we map the normalised root-mean-square (r.m.s.) velocity fluctuations,
$u_{\mathit{rms}}'/\langle u \rangle _{t,\mathit{\boldsymbol{x}}}$
, where
$\langle u \rangle _{t,\mathit{\boldsymbol{x}}}$
is the mean velocity magnitude averaged over time
$t$
and space
$\boldsymbol{x}$
, and the spatial field of r.m.s. velocity fluctuations is
$u_{\mathit{rms}}'(\boldsymbol{x})=\sqrt {\langle (u(\boldsymbol{x},t)-\langle u(\boldsymbol{x})\rangle _t)^2+ (v(\boldsymbol{x},t)-\langle v(\boldsymbol{x})\rangle _t)^2\rangle _t}$
. The fluctuations corresponding to the plane imaged in figure 3(b) are shown in figure 4(a). At
$\textit{Wi} = 1.5$
, the flow is steady and thus no appreciable fluctuations are seen. Above
$\textit{Wi}=\textit{Wi}_c\approx 3$
, sustained flow fluctuations larger than the PIV noise threshold (
$u'/\langle |\boldsymbol{u}| \rangle _t \gt 0.2$
) start to appear. At
$\textit{Wi} = 4.5$
, velocity fluctuations are largest in the bulk void space of the unit cell, reflecting the cross-flow events that intermittently wash away eddies. Figure 4(a) also presents kymographs (showing space and time on the vertical and horizontal axes, respectively) of the normalised instantaneous velocity fluctuations for a selected column of pixels located in a region upstream of the inter-grain contacts over a duration
$\tau _{PV}=100$
. The kymograph for
$\textit{Wi} = 1.5$
shows no fluctuations over space and time because the flow is steady, whereas the kymograph for
$\textit{Wi}=4.5$
shows velocity fluctuations persisting over both space and time. The dark purple bands in the
$\textit{Wi} = 4.5$
kymograph correspond to the eddy-dominated regions, and eddy-washing events are observed as a thinning of these dark purple bands over time. Figure 4(b) presents maps of the normalised r.m.s. fluctuations for additional
$z$
-planes (panels i–viii) spanning half of the characteristic unit cell volume for
$\textit{Wi} = 4.5$
. The largest velocity fluctuations are located in the bulk void space and the
$z$
-plane in panel (viii). These
$z$
-planes correspond to regions of the unit cell containing contacts that are directly accessible to the impinging flow. Note that the neck-like contacts in figure 4(b-viii) are rotated
$90^\circ$
with respect to the contacts in the
$z$
-plane shown in figure 4(a); although these contacts are rotated with respect to the imaging plane, they generate equivalent stagnation points to those in the
$z$
-plane presented in figure 4(a). Altogether, our flow visualisation demonstrates that the additional stagnation points generated by inter-grain contacts causes the elastic instability to develop locally near the contacts compared with other regions in the pore volume.
Characterisation of eddy development and spectral analysis of velocity fluctuations in the SC packing. (a) Quantification of eddy regimes: probability density plot of eddy length normalised by the maximum eddy length. For
$\textit{Wi} \lt 2$
, no eddies are present and the flow is steady. For
$2 \lt \textit{Wi} \lt 3$
, the average eddy length increases with
$\textit{Wi}$
and remains fairly steady over time. Above
$\textit{Wi} = 3$
, eddy lengths can reach their maximum, but are also highly fluctuating as reflected by the large spread in the length distributions. Eddy length distributions at each
$\textit{Wi}$
correspond to measurements sampled every 10 s over 5 minutes of imaging. (b) Temporal power spectral densities
$E_v(f)$
for the transverse velocity fluctuations
$v'$
in the SC packing. We determine the temporal power spectrum at each PIV grid location in the field of view, then average the temporal spectra over all positions to reduce noise. The onset of steep power spectral decay coincides with the onset of instability
$\textit{Wi}_c \approx 3$
. Power laws are fitted to the regions of the spectra with the steepest decay. At
$\textit{Wi} = 4.5$
, the decay exponent is
$\alpha = 1.9$
. The grey vertical line marks the inverse of the longest polymer relaxation time,
$1/\lambda$
(
$\lambda = 0.4$
s).
We further characterise this flow instability by measuring the normalised eddy length
$L_{\mathit{eddy}}/L_{\mathit{max}}$
every 10 s over 5 minutes of imaging for each
$\textit{Wi}$
tested. Briefly, we obtain 5 min. movies at 30 frames per second at each
$\textit{Wi}$
and generate pathline movies in ImageJ using a moving average of 30 frames. We then manually measure eddy lengths
$L_{\mathit{eddy}}$
in ImageJ with the origin of the eddy located at the upstream face of the grain contact and sample eddy lengths every 10 s. We normalise the measured lengths by the maximum possible eddy length
$L_{\mathit{max}}$
, which is equal to the streamwise distance between two consecutive grain contacts. We find three distinct flow regimes, as shown by the probability density plot in figure 5(a). For
$\textit{Wi} \lt 2$
, the flow is steady and eddy-free. At larger
$2 \lt \textit{Wi} \lt 3$
, eddies begin to form on the upstream faces of sphere contacts and grow increasingly longer with
$\textit{Wi}$
(supplementary movie 1). In this regime, the eddies remain fairly steady in length, reflected by a narrower spread of the sampled eddy lengths. Above the threshold
$\textit{Wi}=\textit{Wi}_c\approx 3$
, the eddies are highly unsteady, with lengths that fluctuate drastically between nearly zero to as large as the distance between adjacent contacts. These fluctuations are also reflected in the power spectral density
$E_v(f)=\lvert \mathit{FFT}(Y)\rvert ^2/ (L F_s )$
of transverse velocity fluctuations
$v'$
of different temporal frequencies
$f$
; here,
$\mathit{FFT}(Y)$
is the fast Fourier transform of the temporal velocity fluctuation signal
$Y$
,
$L$
is the length of the signal and
$F_s$
is the signal sampling rate. As shown by the green to yellow curves in figure 5(b), above the threshold
$\textit{Wi}_c \approx 3$
, the spectra begin to exhibit a power law decay
$E_v(f)\sim f^{-\alpha }$
at frequencies smaller than
$1/\lambda$
. The total power contained in these fluctuations, given by the area under the curve of
$E_v(f)$
, increases with
$\textit{Wi}$
, with most of the energy contained in frequencies smaller than
$1/\lambda$
: the flow fluctuations occur over a broad range of time scales primarily longer than the longest polymer relaxation time. The characteristic power law exponent
$\alpha \approx 2$
is consistent with previous measurements in model 3-D porous media at similar
$\textit{Wi}$
(Browne & Datta Reference Browne and Datta2021; Carlson et al. Reference Carlson, Toda-Peters, Shen and Haward2022). This exponent is lower than the value of
$\alpha \approx 3.5$
that is typically referenced as a spectral signature of purely ‘elastic turbulence’ (Groisman & Steinberg Reference Groisman and Steinberg2000; Fouxon & Lebedev Reference Fouxon and Lebedev2003; Groisman & Steinberg Reference Groisman and Steinberg2004; Steinberg Reference Steinberg2019; Varshney & Steinberg Reference Varshney and Steinberg2019; van Buel & Stark Reference van Buel and Stark2022), suggesting that our experiments may probe a transitional unsteady flow state en route to fully developed elastic turbulence.
Elastic instability in the BC packing: chaotic, ‘wobbling’ pathlines. (a) Schematics contrasting an idealised BC cuboid unit cell (top) containing pure point contacts with the consolidated BC cuboid unit cell (bottom) used in this work. The small neck-like region outlined in red between spheres marks an example of an inter-grain contact. Red
$\times$
s correspond to stagnation points generated at sphere surfaces. Yellow
$\times$
s correspond to stagnation points arising at the contacts. Pink
$\times$
s also mark contact-generated stagnation points, but for contacts that are rotated 90
$^\circ$
with respect to the
$z$
-plane shown in the pathline images. (b) Pathline images corresponding to the
$z$
-plane containing the inter-grain contacts (circled in red). (i) Steady flow field at
$\textit{Wi} = 1.5$
. Pathlines are generated by averaging frames over 20 s. Scale bar is 200
$\unicode {x03BC}$
m. (ii) Unsteady flow field at
$\textit{Wi} = 3.5$
for three different time points. The instability is characterised by crossing pathlines; the yellow box highlights a region of pronounced pathline wobbling upstream of a contact. Pathlines are generated by averaging frames over 2 s. (c) Pathline images in the
$z$
-plane located 1/2 sphere diameter above the
$z$
-plane in panel (b). (i) Steady flow field at
$\textit{Wi} = 1.5$
. (ii) Unsteady flow field at
$\textit{Wi} = 3.5$
for three different time points. The pore space geometry in this
$z$
-plane is equivalent to translating the field-of-view for the
$z$
-plane in panel (a) due to the symmetry of the BC unit cell.
3.1.2. Body-centred cuboid packing: wobbling pathlines
We next consider flow of the same polymer solution through the BC packing. This medium also has neck-like contacts between adjacent grains (outlined in red in figure 6
a). As we detail later, by generating stagnation points (yellow and pink
$\times$
s in figure 6
a), these contacts shape the unsteady flow behaviour – just as in the SC case. However, the BC packing contains another distinct set of stagnation points: those at the polar upstream/downstream grain surfaces (red
$\times$
s in figure 6
a), as in the case of Stokes flow around a single sphere. As we show below, this change in geometry changes features of how the elastic instability manifests in BC packings.
Figure 6(b) shows pathline images of the flow in a
$z$
-plane containing inter-grain contacts. At low
$\textit{Wi} = 1.5$
, the flow is steady with steady, uncrossing pathlines (panel i). As
$\textit{Wi}$
increases above a threshold value, the elastic instability develops upstream of the contacts (panel ii) – just as in the SC case – rather than at the polar upstream/downstream stagnation points. However, this unsteady flow is different from that in the SC case in two ways. First, the threshold
$\textit{Wi}_c\approx 2$
is slightly smaller than the threshold
$\textit{Wi}_c\approx 3$
in the SC packing. Second, in the BC packing, the elastic instability does not manifest as unsteady eddies, but rather pronounced pathline ‘wobbling’. An example is shown by the crossing pathlines upstream of the contacts shown in the boxed region in figure 6(b-ii) for
$\textit{Wi}=3.5$
. Supplementary movie 4 shows this transition from steady to unsteady flow, with increasingly chaotic pathline wobbling with increasing
$\textit{Wi}$
. The instability is initially generated upstream of the inter-grain contacts before fluctuations start to persist throughout the entire pore space at
$\textit{Wi} \approx 4$
.
Figure 6(c) shows pathline images for the
$z$
-plane located half a grain diameter above that shown in figure 6(b). We observe similar wobbling behaviour upstream of the contacts, which reflects the fact that these two
$z$
-planes are symmetry planes due to the periodicity of the BC packing. Supplementary movies 5 and 6 provide additional pathline imaging for steady (
$\textit{Wi}=1.5$
) and unsteady (
$\textit{Wi} = 3.5$
) flow states across nine intermediate
$z$
-planes that span half of the characteristic unit cell volume. At sufficiently high
$\textit{Wi}$
above the onset of the instability (
$\textit{Wi}_c\approx 2$
), the flow throughout the entire pore volume is unsteady, as shown in supplementary movie 6.
Spatiotemporal features of the BC wobbling instability. (a) Left: Normalised r.m.s. velocity fluctuation colourmaps for steady (
$\textit{Wi}=1.5$
, top left) and unsteady (
$\textit{Wi}=3.5$
, bottom left) flow states. Right: Kymographs of the normalised velocity fluctuations corresponding to the column of pixels boxed in colourmaps to the left. For
$\textit{Wi} = 1.5$
, each column in the kymograph corresponds to a chunked average over 15 frames to reduce noise. No averaging is applied to the
$\textit{Wi} = 3.5$
kymograph. Time
$t$
is reported as the number of injected pore volumes,
$\tau _{PV}=t/(\mathit{PV}/Q)$
, where one pore volume (PV) is defined as the void volume of one unit cell of the packing. (b) Normalised r.m.s. velocity fluctuation maps for additional
$z$
-planes at
$\textit{Wi}=3.5$
spanning half of the unit cell. Sequential
$z$
-planes are spaced approximately
$50 \: \unicode {x03BC}$
m apart.
Temporal power spectral densities
$E_v(f)$
for the transverse velocity fluctuations
$v'$
in the BC packing at various
$\textit{Wi}$
. The onset of power spectral decay coincides with the onset of instability,
$\textit{Wi}_c \approx 2$
. Power laws are fitted to the regions of the spectra with the steepest decay. The grey vertical line marks the inverse of the longest polymer relaxation time (
$\lambda = 0.4$
s).
We again map the r.m.s. velocity fluctuations normalised by the mean flow speed,
$u_{\mathit{rms}}'/\langle u \rangle _{t,\mathit{\boldsymbol{x}}}$
, as shown in figure 7(a). At
$\textit{Wi} = 1.5$
, the flow is steady with no noticeable fluctuations. At
$\textit{Wi} = 3.5$
above the onset of instability, pronounced velocity fluctuations occur upstream of the inter-grain contact, concomitant with the unsteady wobbling behaviour seen in the pathline imaging (figures 6
b and 6
c). Figure 7(a) also presents kymographs of the normalised instantaneous velocity fluctuations
$u'/\langle u \rangle _{t,\mathit{\boldsymbol{x}}}$
over 50 injected pore volumes for a column of pixels located upstream of the contact area. No spatiotemporal fluctuations occur at
$\textit{Wi} = 1.5$
, contrasted by bursts of large fluctuations at
$\textit{Wi} = 3.5$
that display no characteristic regularity over time; consequently, the flow does not exhibit a periodic state during the transition to instability. Figure 7(b) provides maps of the normalised r.m.s. velocity fluctuations for unsteady flow at
$\textit{Wi}=3.5$
across additional
$z$
-planes spanning half of the characteristic unit cell volume. Fluctuations in the planes containing the contacts (figures 7
a and 7b-viii) are largest upstream of the contacts, while fluctuations in the intermediate planes (figures 7
b-ii–vii) are less localised throughout the pore space, perhaps due to the more uniform curvature of the geometry.
Figure 8 presents the temporal power spectral densities
$E_v(f)$
of the transverse velocity fluctuations
$v'$
in the BC packing at various
$\textit{Wi}$
. As in the SC case, the fluctuations are persistent: most of the kinetic energy in the velocity fluctuations is contained in frequencies smaller than
$1/\lambda$
. Also, the total power contained in these fluctuations increases with
$\textit{Wi}$
. Finally, the spectra similarly exhibit a broad power law decay
$E_v(f)\sim f^{-\alpha }$
, with
$\alpha \approx 2$
, above the threshold
$\textit{Wi}_c = 2$
– just as in the SC packings and consistent with previous measurements at similar
$\textit{Wi}$
in 3-D porous media (Browne & Datta Reference Browne and Datta2021; Carlson et al. Reference Carlson, Toda-Peters, Shen and Haward2022). As in the SC case, this exponent suggests that our experiments may probe a transitional state that arises en route to fully developed purely elastic turbulence.
3.2. Local Pakdel–McKinley criterion rationalises the role of inter-grain contacts in the onset of the elastic instability
Our experiments have demonstrated that, in both the SC and BC geometries, the spatially extended contacts formed between adjacent grains generate stagnation points that lead to the development of the elastic instability. For the BC packing, which also has stagnation points at the upstream/downstream poles of grains, the points associated with grain contacts instead generate the instability. As we show next, these findings can be rationalised within the conventional theoretical framework that has been used to understand the onset of elastic instabilities in other, simpler geometries.
In particular, to account for the role of both streamline curvature and fluid elasticity in destabilising the flow, McKinley et al. (Reference McKinley, Pakdel and Öztekin1996) proposed a dimensionless parameter, the
$M$
parameter, that aims to describe the onset of elastic instabilities for different solution rheology and flow geometries. The
$M$
parameter is defined as
$M=\sqrt {({\lambda U}/{\mathcal{R})} ({N_1}/{\eta _0 \dot {\gamma }})}$
, where
$U$
is the characteristic flow speed and
$\mathcal{R}$
is the streamline radius of curvature; it thus requires both sufficient streamline curvature and elastic normal stresses for the instability to occur. Critical values of
$M_c \approx 1-10$
have been reported for a range of model geometries (McKinley et al. Reference McKinley, Pakdel and Öztekin1996; Pakdel & McKinley Reference Pakdel and McKinley1996), although a universal transition value of
$M$
has not been found experimentally (Casanellas et al. Reference Casanellas, Alves, Poole, Lerouge and Lindner2016). To assess whether this parameter can help describe the onset of the elastic flow instability in our experiments, we use the flow fields measured in both packings to directly compute
$M=\sqrt {(N_1(\dot {\gamma })/\eta _0 \dot {\gamma }) \: (\lambda |\mathit{u}|\kappa )}$
, where
$\dot {\gamma }$
is the magnitude of the rate-of-strain tensor
$|\boldsymbol{D}|=\sqrt {2\boldsymbol{D}:\boldsymbol{D}}$
,
$\eta _0$
is the fluid zero-shear viscosity,
$|\mathit{u}|$
is the local velocity magnitude and
$\kappa$
is the local streamline curvature, which we determine as a 2-D local streamline curvature. Under laminar flow conditions, the symmetry of the ordered sphere packings requires that the flow fields in most of the interrogated
$z$
-planes are 2-D; thus, we use the local in-plane streamline curvature determined from experimental 2-D velocity fields to evaluate the local
$M$
parameter.
Colourmaps of the local
$M$
parameter normalised by the maximum
$M$
value for the (a) SC (
$\textit{Wi} = 2.5$
) and (b) BC (
$\textit{Wi} = 1.5$
) packings at flow conditions below the onset of instability. The largest relative
$M$
values for both sphere packings are located in the pore space near the inter-grain contacts, suggesting that the contacts generate sufficient streamline curvature and elastic stress for the instability to develop. The local
$M$
fields are temporally averaged over 2 minutes of imaging to reduce noise. Maps are presented for nine
$z$
-planes (from left to right) spanning half of the characteristic unit cell and sequential
$z$
-planes are spaced approximately
$50 \: \unicode {x03BC}$
m apart. Example contacts are marked by black circles. The maximum values of M at these flow conditions are
$M_{\mathit{max}}= 1.7$
and 0.82 for the SC and BC packings, respectively.
Figure 9 shows the resulting colourmaps of
$M$
normalised by its maximum value,
$M(\boldsymbol{x})/\mathit{max}(M)$
, taken at
$\textit{Wi}$
values just below the instability threshold for nine different
$z$
-planes spanning half of the characteristic unit cell volume. For the SC packing (figure 9
a), the regions of highest
$M$
are located in the
$z$
-planes nearest the inter-grain contacts. For the BC packing (figure 9
b), the regions of highest local
$M$
occur as lobes on either side of the contacts. The origin of the lobes in the BC packing arises from the dominant contribution of large streamline curvature in fluid streamlines traversing around the sides of the grain contacts. Indeed, in both geometries, the highest streamline curvature under steady flow conditions occurs locally at the contacts, which have a maximum curvature
$\kappa \sim 1/R_{\mathit{neck}}$
. These measurements demonstrate that the flow around inter-grain contacts – rather than around the upstream/downstream poles traditionally considered for isolated spheres – has the largest streamline curvature and, therefore, these contacts control the onset of the elastic instability as described theoretically by McKinley et al. (Reference McKinley, Pakdel and Öztekin1996).
Common transition to elastic instability. (a) Fraction of time that a pore is unsteady (
$F_t$
) plotted against the macroscopic
$\textit{Wi}$
rescaled by the pore-specific critical Weissenberg number
$\textit{Wi}_{c,\mathit{pore}}$
for randomly selected pores across each porous medium, where the instability threshold criterion is
$u'/|\boldsymbol{u}| \gt 0.2$
. The transition to instability collapses for the ordered packings and disordered sphere packing. The solid black line is consistent with power law scaling
$F_t\,\sim (\textit{Wi}/\textit{Wi}_c-1)^\alpha$
with
$\alpha = 0.4 \pm 0.1$
, although more measurements are needed to determine the exact nature of this transition. (b) Distributions of
$\textit{Wi}_{c,\mathit{pore}}$
for randomly selected pores across each medium. The distributions are much narrower for the SC and BC geometries compared with the disordered packing, consistent with the hypothesis of porous individualism. Data corresponding to the disordered packing are reproduced from Browne & Datta (Reference Browne and Datta2021).
3.2.1. Transition to instability collapses across different porous medium geometries
Given that the framework of McKinley et al. (Reference McKinley, Pakdel and Öztekin1996) provides a common way to describe the onset of the elastic instability across different geometries in our experiments, we next examine whether there may be any commonalities in how the flow transitions from steady to unsteady states. We characterise the extent of this transition by determining the fraction of time that pores are unsteady
$F_t$
, where values of
$0\lt F_t\lt 1$
correspond to transitory flow states with pores switching between steady and unsteady states over time. We image 10 randomly selected pores in each medium and apply an instability threshold of
$u'/\langle |\boldsymbol{u}| \rangle _t \gt 0.2$
for at least 10 % of pixels in the frame to determine the instantaneous stability of each pore. This stability threshold was chosen to ensure that fluctuations calculated from the flow field are above the noise threshold of particle image velocimetry. The instantaneous velocity fluctuation fields are calculated from
$u'(\boldsymbol{x},t)=\sqrt {((u(\boldsymbol{x},t)-\langle u(\boldsymbol{x}\rangle _t)^2+(v(\boldsymbol{x},t)-\langle v(\boldsymbol{x}\rangle _t)^2}$
. Figure 10(a) plots
$F_t$
against the rescaled Weissenberg number for each pore,
$\textit{Wi}/\textit{Wi}_{c,\mathit{pore}}-1$
, where
$\textit{Wi}_{c,\mathit{pore}}$
is the threshold
$\textit{Wi}_c$
for individual pores as determined by the instability criterion detailed previously. We also show data for the flow of an identical HPAM solution through a disordered random sphere packing with grain sizes of
$D_{\!p} = 300-355 \: \unicode {x03BC}$
m and porosity
$\phi \approx 0.41$
, reproduced from Browne & Datta (Reference Browne and Datta2021). The transition plots appear to collapse, albeit with some spread among the data, both across pores within the same geometry and across the three geometries overall, over a narrow band of
$\textit{Wi}/\textit{Wi}_c$
– suggesting a common pathway to instability despite the different packing types. However, more measurements are needed near
$\textit{Wi}_c$
to determine precise scaling exponents to characterise the exact nature of the transition and the resulting elastic instability.
Previously, Browne & Datta (Reference Browne and Datta2021) proposed the concept of ‘porous individualism’: pores of different geometries in a disordered porous medium can exhibit different flow states for the same imposed flow rate. This feature is reflected by a fairly broad distribution of threshold
$\textit{Wi}_{c,\mathit{pore}}$
for individual pores sampled in a disordered packing, as shown in figure 10(b) (reproduced from Browne & Datta (Reference Browne and Datta2021)). By contrast, unlike the random disordered packing, all pores in the ordered packings studied here have nearly identical geometries. Indeed, we measure much narrower distributions of
$\textit{Wi}_{c,\mathit{pore}}$
for pores sampled across each medium, as shown in figure 10(b). We hypothesise that the width of the distributions results from polymer stress history effects, where variations in the polymeric stress field exist from pore-to-pore due to fluid memory. Additionally, small variations in
$\textit{Wi}_c$
could arise from experimental imperfections resulting from the fabrication process. Further, we observe that the SC and disordered packing
$\textit{Wi}_{c,\mathit{pore}}$
distributions overlap, while the BC distribution is located at overall lower
$\textit{Wi}_{c,\mathit{pore}}$
values. It is still unclear what sets
$\textit{Wi}_{c,\mathit{pore}}$
for each geometry and investigating this point further will be an important direction for future work. However, we conjecture that this difference in onsets for the BC and SC geometries results primarily from a difference in the net stagnation point densities: the BC packing has a higher net overall density due to the presence of both inter-grain contact and upstream/downstream polar stagnation points.
3.3. Flow topology reflects extensionally dominated flow due to stagnation points
Our analysis thus far has assumed that stagnation points – which are dominated by locally extensile flow – are generated at inter-grain contacts and these points generate the elastic instability at sufficiently large
$\textit{Wi}$
. To directly confirm this assumption, we use the measured flow fields to determine the local flow type parameter
$\xi =(|\boldsymbol{D}|-|\boldsymbol{\varOmega }|)/(|\boldsymbol{D}|+|\boldsymbol{\varOmega }|)$
, where
$|\boldsymbol{\varOmega }|=\sqrt {2\boldsymbol{\varOmega }:\boldsymbol{\varOmega }}$
is the magnitude of the vorticity tensor (Astarita Reference Astarita1979; Hamlington, Schumacher & Dahm Reference Hamlington, Schumacher and Dahm2008). The flow type parameter takes on values from −1 to 1:
$\xi =-1$
corresponds to solid-body rotation,
$\xi =0$
corresponds to simple shear and
$\xi = +1$
corresponds to pure extension. Figure 11 presents maps of
$\xi$
for nine
$z$
-planes spanning half of the characteristic unit cell volume at steady flow conditions (
$\textit{Wi} = 1.5$
) in both packings, and figures 12(a), 12(b) present full distributions of the measured values of
$\xi$
that quantify the relative occurrences of flow type within the pore volume of each packing.
Local flow type parameter (
$\xi$
) maps for the (a) SC and (b) BC packings at steady flow conditions (
$\textit{Wi} = 1.5$
). Extensional flow is dominant near the inter-grain contacts, and regions of the flow field with expansion and contraction of streamlines. Colourmaps are presented for nine
$z$
-planes spanning half of the characteristic unit cell volume, where sequential
$z$
-planes are spaced approximately
$50 \: \unicode {x03BC}$
m apart. Examples of contacts and contact-generated stagnation points are marked by black circles and yellow
$\times$
s, respectively. The local
$\xi$
fields are temporally averaged over 2 minutes of imaging to reduce noise. Scale bar is
$200 \: \unicode {x03BC}$
m.
Quantification of the extensional nature of the flow in the ordered packings. (a,b) Probability density functions (p.d.f.s) of the flow type parameter,
$\xi$
, at
$\textit{Wi} = 1.5$
in the SC and BC packings. (c,d) Complementary cumulative distribution functions (1-CDFs) of the Hencky strain,
$\epsilon _{\mathit{Hencky}}$
, accumulated over 1 relaxation time across different
$\textit{Wi}$
in the SC and BC packings. Appreciable values of the Hencky strain are accumulated in each geometry above the respective critical
$\textit{Wi}$
, with the largest strains accumulating in the BC packing. Each distribution comprises 75 values.
In the SC packing (figure 11
a), extensional flow dominates upstream and downstream of the inter-grain contacts, as expected. Extensional flow also exists in the middle of the bulk void space due to the contraction and expansion of the streamlines as the fluid navigates between successive pore bodies. These regions of local extension, arising from the stagnation points, are reflected by the skewed tail towards extensional flow shown in the measured distribution of
$\xi$
(figure 12
a). Similarly, in the BC packing (figure 11
b), extensional flow is found upstream and downstream of both the contacts and the upstream/downstream poles of the grains, while flow in the intermediate planes appears to be more shear-dominated due to the increased confinement. As a result, the flow topology distribution in the BC packing is more strongly skewed towards extensional flow (
$\xi = 1$
) with a peak also occurring near
$\xi =0$
(figure 12
b). Altogether, this in situ visualisation directly confirms that inter-grain contacts generate extension-dominated stagnation points.
While the flow-type parameter can be used to identify regions of different flow topology, it does not provide information on the strength of each flow type. Instead, to quantify the strength of extensional effects, we calculate the Hencky strain accumulated by a polymer over the characteristic relaxation time using the local extension rate from the 2-D velocity field,
$\dot {\epsilon }=|{\partial u}/{\partial x}|+|{\partial v}/{\partial y}|$
. The Hencky strain describes the net extensional deformation of a polymer chain due to an extensional velocity gradient over a certain time period; for sufficiently large Hencky strains, the polymers can undergo a coil–stretch transition (De Gennes Reference De Gennes1974), whereby the polymers become elongated and large extensional stresses develop (Anna et al. Reference Anna, McKinley, Nguyen, Sridhar, Muller, Huang and James2001). We select a starting location in the pore space and integrate forwards in time along the Lagrangian trajectory to determine the accumulated Hencky strain,
$\epsilon _{\mathit{Hencky}}\approx \sum _{t=0}^{\lambda }\dot {\epsilon }\Delta t$
. We repeat this for five random positions at 15 different time points to generate distributions of
$\epsilon _{\mathit{Hencky}}$
at different
$\textit{Wi}$
. Figures 12(c) and 12(d) present the resulting complementary cumulative distribution functions for the SC and BC packings. In both geometries,
$\epsilon _{\mathit{Hencky}}$
increases with
$\textit{Wi}$
, reaching larger values in the latter compared with the former at matching
$\textit{Wi}$
– reflecting the greater incidence of extension-dominated stagnation points in the BC packing. Moreover, for
$\textit{Wi}\gt \textit{Wi}_c$
in both packings,
$\epsilon _{\mathit{Hencky}}$
predominantly exceeds 2 – the value at which the extensional viscosity of dilute polymer solutions is known to drastically increase as a result of polymer stretching (Anna et al. Reference Anna, McKinley, Nguyen, Sridhar, Muller, Huang and James2001; McKinley & Sridhar Reference McKinley and Sridhar2002). In contrast, for the flow of the same polymer solution in a disordered random sphere packing at similar values of
$\textit{Wi}$
, the magnitudes of the accumulated Hencky strains are smaller than 1 (Browne & Datta Reference Browne and Datta2021). Therefore, we expect flow in the ordered packings to be more extensionally dominant compared with more shear-dominated flow in a random sphere packing. As we show next, this difference in the flow type and resulting difference in polymer extension has marked implications for the macroscopic flow resistance through a porous medium.
3.4. Flow thickening is generated by the elastic instability, primarily due to extensional effects
Previous studies of HPAM solution flow resistance in porous rock cores (Seright et al. Reference Seright, Fan, Wavrik and Balaban2010; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland, Hawkes and Leeper2015; Mitchell et al. Reference Mitchell, Lyons, Howe and Clarke2016) and disordered 3-D sphere packings (Browne & Datta Reference Browne and Datta2021) reported anomalous flow thickening – marked by a steep increase of the apparent viscosity
$\eta _{\mathit{app}}$
above the bulk shear viscosity
$\eta _I$
– for
$\textit{Wi}$
increasing above a threshold value. Browne & Datta (Reference Browne and Datta2021) demonstrated that this behaviour reflects the onset of an elastic instability and established that the increase in measured macroscopic flow resistance arises from the added viscous dissipation from fluctuations in the unsteady flow field. However, their study was for a disordered medium, in which the flow is more shear-dominated, whereas our measurements here reveal that flow in ordered 3-D packings is more extensionally dominant due to stagnation points generated at inter-grain contacts. What are the implications of this difference in flow type for the macroscopic flow resistance?
Macroscopic flow resistance increases above a threshold
$\textit{Wi}$
for the same polymer solution in different 3-D sphere packing geometries. The apparent viscosity
$\eta _{\mathit{app}}$
determined using Darcy’s Law is normalised by the viscosity measured from bulk shear rheology
$\eta _I$
. The dotted line at
$\eta _{\mathit{app}}/\eta _I = 1$
is the Darcy flow resistance predicted by the bulk shear rheology. Solid curves correspond to the power balance model incorporating an added viscous dissipation generated by unsteady flow fluctuations. Shaded bands reflect the uncertainty in the power balance model. Data for the disordered packing are reproduced from Browne & Datta (Reference Browne and Datta2021). The power balance captures the onset of flow thickening in the ordered packings, but increasingly underpredicts the flow resistance with increasing
$\textit{Wi}$
. Discrepancies between the power balance model and the measured flow resistance likely arises from extension-induced polymer stresses generated by contact-generated stagnation points in ordered packings, as indicated by the double-headed arrows. In some cases, error bars may be smaller than symbol size.
To address this question, we measure the pressure drop across each medium as a function of
$\textit{Wi}$
and determine the corresponding apparent viscosity
$\eta _{\mathit{app}}$
using Darcy’s law. We define a non-dimensional measure of flow resistance by normalising
$\eta _{\mathit{app}}$
by the bulk shear viscosity
$\eta _I$
evaluated at the corresponding interstitial shear rate. As shown in figure 13, both the SC and BC packings exhibit flow thickening above a threshold
$\textit{Wi}=\textit{Wi}_c$
, with
$\textit{Wi}_c\approx 3$
and
$\textit{Wi}_c\approx 2$
for SC and BC, respectively – coincident with the onset of the elastic instability (figure 10
b), as was found previously for disordered media as well. Interestingly, we do not observe a downturn in the flow resistance at high
$\textit{Wi}$
as previously observed in disordered sphere packings (Browne & Datta Reference Browne and Datta2021); given the limited range of
$\textit{Wi}$
accessible in our experiments, it will be useful in future work to explore the flow resistance behaviour at high
$\textit{Wi}$
where reduced flow-thickening or even drag reduction behaviour may be observed (Savins Reference Savins1969; Müller & Sáez Reference Müller and Sáez1999; Skauge et al. Reference Skauge, Zamani, Gausdal Jacobsen, Shaker Shiran, Al-Shakry and Skauge2018; Varshney & Steinberg Reference Varshney and Steinberg2018; Kumar et al. Reference Kumar, Varshney, Li and Steinberg2022b
).
Our in situ flow visualisation also provides a way to quantify the extent to which this flow thickening is due to the added fluid viscous dissipation from fluctuations generated by the instability, following the framework of Browne & Datta (Reference Browne and Datta2021). Specifically, we directly compute the time-averaged rate of added viscous dissipation per unit volume due to pore-scale fluctuations from the flow field measured in single unit cells:
$\langle \chi \rangle _{t,V}=\eta _I\langle \boldsymbol {D}':\boldsymbol {D}'\rangle _{t,V}$
, where
$\boldsymbol{D}'$
is the rate-of-strain tensor associated with velocity fluctuations from the mean and we average over nine
$z$
-planes spanning half of the characteristic unit cell volume. We fit a power law model to describe the growth of the added viscous dissipation with
$\textit{Wi}$
:
$\langle \chi \rangle _{t,V}=\alpha (\textit{Wi}-\textit{Wi}_c)^\beta$
, where the corresponding parameters are
$\alpha _{\mathit{SC}}=7.4\pm 0.5 \ \mathrm{Wm}^{-3}$
,
$\alpha _{\mathit{BC}}=2.3\pm 0.2 \ \mathrm{Wm}^{-3}$
,
$\beta _{\mathit{SC}}=1.5\pm 0.1$
and
$\beta _{\mathit{BC}}=1.3\pm 0.1$
. As detailed by Browne & Datta (Reference Browne and Datta2021), balancing the rate of work done by the fluid pressure with the rate of viscous energy dissipation per unit volume yields an expression for how this added dissipation influences macroscopic flow resistance:
\begin{equation} \eta _{\mathit{app}}/\eta _I = 1 +k \langle \chi \rangle _{t,V}/(Q/A_c)^2 \eta _I. \end{equation}
The rightmost term in (3.1) quantifies the resistance from the added viscous dissipation generated by the elastic instability and functionally depends on
$\textit{Wi}$
as
${\sim} ({(\textit{Wi}-\textit{Wi}_c)^\beta }/{\textit{Wi}^{(n+1)/(n_{N_1}-n)})}$
. This analysis assumes a power law fluid model that captures shear-thinning but neglects polymer elasticity. Therefore, it quantifies the additional power dissipated in the solvent due to flow fluctuations, but neglects the contribution of the polymeric stresses themselves, which we expect to be increasingly dominant in the ordered packings due to the greater prevalence of extensional flow.
This expectation is borne out by our data. Evaluating (3.1) using the measured flow fields produces the solid curves shown in figure 13. For both of the ordered packings, (3.1) captures the onset of flow thickening, but increasingly underpredicts the flow resistance at higher
$\textit{Wi}$
– consistent with our expectation that extension-induced polymeric stresses play a more dominant role in these ordered sphere packings. Indeed, our local measurements of Hencky strain and flow type parameter highlight differences between extensionally dominated flows in ordered packings and more shear-dominated flows in random, disordered packings, ultimately hinting at different mechanisms of increased flow resistance in ordered and disordered geometries. Moreover, the discrepancy in the power balance from the measured flow resistance is even larger for the BC packing, for which we found that extensional flow is even greater (figures 11 and 12). Altogether, these data suggest that extensional viscosity arising from stagnation point-induced polymer extension has a larger contribution to the flow resistance in the SC and BC geometries, particularly with increasing
$\textit{Wi}$
. Incorporating a spatially varying polymer stress distribution into the model may improve the predictive power of (3.1) and will be a useful direction for future work; however, in situ measurements of local polymer stress fields in porous media remain challenging experimentally.
4. Conclusions
Our work has revealed that the spatially extended contacts between grains in a 3-D medium, which are typically idealised as discrete points and therefore overlooked, play a pivotal role in shaping the unsteady flow of a polymer solution through the pore space. They generate stagnation points that drive the onset of an elastic instability in both SC and BC sphere packings, contrasting the hypothesis that polar upstream and downstream stagnation points at sphere surfaces are the primary generator of elastic instabilities in flow past spheres. This finding suggests that not only is the effective stagnation point density within a particular geometry important, as suggested by Haward et al. (Reference Haward, Hopkins and Shen2021) for 2-D pillar arrays, but that different classes of stagnation points may arise in flows through 3-D porous media dependent on the porous medium geometry. Indeed, stagnation points arising at inter-grain contacts are inherently inaccessible in 2-D model geometries, which many lab studies employ, yet are intrinsic to many natural porous media.
A simple rearrangement of the sphere packing type changes the pore space geometry confining the fluid, altering the distribution and density of stagnation points that contribute to the development of the elastic instability as well as the underlying flow topology. As a result, we expect that the features of the instability will be different for flow of the same polymer solution in the two different sphere packing geometries. Our in situ visualisation confirmed this expectation: the instability generates fluctuating 3-D eddies in the SC packing versus wobbling streamlines in the BC packing. Despite this marked difference in the pore-scale features of the unsteady flow between medium geometries, we found that there are common features of the transition pathway to instability across both the ordered and disordered porous media – suggesting that the route to elastic instability may be common across different pore structures. Thus, more careful examination of the precise nature of the transition to instability and the instability itself will be a valuable direction for future work.
Haward & Odell (Reference Haward and Odell2003) and Odell & Haward (Reference Odell and Haward2006) examined flow of dilute polymer solutions through model simple cubic and body-centred sphere packings, reporting an increase in macroscopic flow resistance that depends on the nature of the polymer, solution properties and packing geometry. However, these and other previous characterisations of flow in 3-D porous media were limited to quantifying bulk flow resistance, as access to in situ flow fields with high spatiotemporal resolution is typically hindered by the opacity of the medium. Our in situ flow visualisation overcame this limitation and revealed the key role of inter-grain contact-generated stagnation points in controlling this flow thickening. In particular, we found that the steep increase in flow resistance coincides with the onset of local velocity fluctuations created by the elastic instability. However, unlike for flow in model pore constriction and expansion arrays (De et al. Reference De, Kuipers, Peters and Padding2017a , Reference De, Kuipers, Peters and Paddingb , Reference De, Kuipers, Peters and Paddingc ) or disordered sphere packings (Browne & Datta Reference Browne and Datta2021; De et al. Reference De, Kuipers, Peters and Padding2017b ), added viscous dissipation arising from unsteady flow fluctuations does not fully capture the flow thickening behaviour in the ordered packings. Rather, our data suggest a dominant role of polymeric stresses generated by the contact-associated stagnation points, which may manifest in complex ways (Chauveteau, Moan & Magueur Reference Chauveteau, Moan and Magueur1984; Evans, Shaqfeh & Frattini Reference Evans, Shaqfeh and Frattini1994; Rothstein & McKinley Reference Rothstein and McKinley2001; Haward & Odell Reference Haward and Odell2003; James, Yip & Currie Reference James, Yip and Currie2012; Zamani et al. Reference Zamani, Bondino, Kaufmann and Skauge2015; Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022). Altogether, our results suggest that disordered porous medium geometries produce shear-dominated flows, consistent with a flow resistance largely governed by viscous dissipation in the solvent (Browne & Datta Reference Browne and Datta2021), whereas ordered porous medium geometries generate more extensionally dominated flow and thus elastic tensile stress-dominated flow resistance. Our current experimental methods lack the capability to quantify in situ polymer stress history; future work focusing on integrating rheological methods and simulations to incorporate stress history and holistically assess multiple mechanisms contributing to flow resistance will be important for fully understanding the mechanisms behind the anomalous flow thickening behaviour of polymer solution flow in porous media.
In this work, we examined ordered porous media to help elucidate the relationship between 3-D pore space geometry and viscoelastic flow behaviour. Although random sphere packings provide a more realistic model of many porous media, it is more challenging to identify the relevant stagnation points in a random packing due to the irregular orientation of grains and pore-to-pore variations in mean flow direction. Investigating local flow behaviour near inter-grain contacts in consolidated, random packings will be a useful direction for future work, as our results suggest that contact area-generated stagnation points have important implications for viscoelastic flow in random porous media. Further, a broad range of porous medium geometries exist in both natural settings and industrial applications beyond the binary classification of ordered versus disordered geometries. Our work focuses on granular-type packings containing consolidations modelled by neck-like contact areas, while other packing morphologies may include loose granular packings or media consisting of continuous solid matrices. We expect these different pore space geometries to fundamentally alter the flow behaviour by controlling the flow topology and stagnation point densities. Indeed, for laminar Stokes flow – in the absence of instabilities altogether – variations in pore space geometry alone can control the onset of chaotic mixing of scalars (Le Borgne, Dentz & Villermaux Reference Le Borgne, Dentz and Villermaux2015; Lester et al. Reference Lester, Dentz and Le Borgne2016; Turuban et al. Reference Turuban, Lester, Le Borgne and Méheust2018, Reference Turuban, Lester, Heyman, Borgne and Méheust2019). Examining polymer solution flow behaviour in controlled model geometries of broader packing morphologies will thus be a useful extension of our work.
By modelling features of naturally occurring porous media in controlled microfluidic systems, our work provides an improved understanding of how pore space geometry influences the onset and features of viscoelastic flow instabilities in 3-D porous media. Our results suggest the tantalising possibility that elastic instabilities could be tailored via design of the pore space geometry for use in targeted applications, such as exploiting eddies to enhance solute mixing (Browne et al. Reference Browne, Shih and Datta2020a ; Chen & Datta Reference Chen and Datta2025) or leveraging wobbling streamlines to promote transverse dispersion of colloids (Scholz et al. Reference Scholz, Wirner, Gomez-Solano and Bechinger2014). Additionally, our results provide intuition for setting operating guidelines where viscoelastic porous media flows are relevant in which flow instabilities can either be prevented or used depending on application objectives.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11458 and https://github.com/emchen6639/lattices.
Acknowledgements
We thank K. Toda-Peters for fabricating the packing devices; D. W. Carlson, A. Clarke and G. H. McKinley for insightful discussions; and the Princeton Imaging & Analysis Center for use of the rheometer.
Funding
We acknowledge support from the Princeton Center for Complex Materials (PCCM), a National Science Foundation (NSF) Materials Research Science and Engineering Center (MRSEC; DMR-2011750, as well as the Camille Dreyfus Teacher-Scholar Program of the Camille and Henry Dreyfus Foundation. S.J.H. and A.Q.S. gratefully acknowledge the support from the Okinawa Institute of Science and Technology Graduate University (OIST), which is funded by a subsidy from the Cabinet Office, Government of Japan. The authors also extend their thanks to the Engineering Section of the Core Facilities at OIST for their valuable assistance. This work was further supported by the Japan Society for the Promotion of Science (JSPS) through Kakenhi grants no. 24K07332 (S.J.H.) and no. 24K00810 (A.Q.S. and S.J.H.).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All data are available in the manuscript and supplementary movies.
Author contributions
Conceptualisation: C.A.B., S.J.H., S.S.D.; funding: S.J.H., A.Q.S., S.S.D.; experiments and data analysis: E.Y.C.; contribution of experiment resources: S.J.H., A.Q.S., S.S.D.; interpretation of results: all authors; manuscript writing: E.Y.C., S.S.D; manuscript editing: all authors.
