Non-Stokesian flow through ordered thin porous media imaged by tomographic-PIV

The 3D flow-fields in a staggered and cubic arrangement of mono-radii cylinders are investigated using tomographic-PIV. The cylinder Reynolds-number is in the range of ≈10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 10$$\end{document} to ≈800\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 800$$\end{document} giving an almost complete overview of the transition region. Two pore-scale effects are discovered. The first, visible in the cubic packing, is a spatially alternating lateral velocity field, which has a significant impact on the pressure drop and transversal dispersion. The second effect, present in the staggered array, is an example of a disturbance propagation effect that takes place in the laminar steady region; this manifests as a peculiar and complex flow-pattern. In accordance with other studies, it is shown that Darcy’s law can, from an engineering point of view be valid far beyond the limit for Stokesian flow.


Introduction
Historically, research on porous media has focused on finding simplified models for macroscopic properties including pressure drop, dispersion and, heat transfer. The development of more advanced measurement and simulation techniques has enabled investigation of the mechanisms leading to these macroscopic properties. Examples of such studies are among others (Comiti and Renaud 1989;Koch and Ladd 1997;Koch and Hill 2001) and (Seguin et al. 1998). For these type of studies the investigation is focused on how the pore-scale dynamics cause discernible macroscopic properties.
As defined by Ziółkowska and Ziółkowski (1988) and Koch and Ladd (1997), there are four distinct flow regimes of interest when studying porous materials, these are distinguished by different flow properties and are defined in terms of the Reynolds number. The Reynolds number defined by Ziolkowska & Ziolkowski is Re z =̄u d z where ū is the mean stream-wise velocity of the liquid within the porous medium, d z is the pore size and is the kinematic viscosity. The different regimes are then approximately: -Stokes flow: Re z << 1 -Laminar flow: 1 < Re z < 10 -Transitional flow: 10 < Re z < 300 -Turbulent flow: Re z > 300 1 3 46 Page 2 of 12 As Ziolkowska & Ziolkowski argues, it is impossible to precisely define the critical Reynolds numbers separating these regions because of the broad variance of reported values from experiments. Therefore the change in flow regime has to be evaluated for each packing separately. The flow regimes affect the macroscopic properties substantially; their characteristics are summarized below.
The Stokes regime where Re z << 1 is characterized by Darcy's law, in equation form p ,i = k ij u j where k ij , in general, is a symmetric 3x3-tensor for 3D-Flows as described by Liakopoulos (1965). For an isotropic porous media the equation reduces to the relation between pressure and flowrate/velocity as.
In other words, the pressure drop increases linearly with flow rate, and as demonstrated by , this law is valid up to Re z ≈ 10 for a packing of spheres. Investigation of the pressure dependence on flow rate has led to Forcheimers postulation of a modified Darcy's law of the form presented in Eq. (2).
If Re z is relatively small, then m = 3 as described by Mei and Auriault (1991). Note that the k 1 factor is vanishingly small for small Reynolds numbers; therefore, Darcy's law is accurate throughout this flow regime. Although Darcy's law is valid throughout the Stokes regime it is not only valid in this regime, for example as Lundstrom et al. 2010;Koch and Ladd 1997) shows it is a good approximation for low Reynolds-number non-Stokesian flows.
The Laminar regime is characterized by steady flow features with inertial effects being non-negligible. Forcheimer's law is still applicable, but the value of m varies with the Reynolds number since the symmetries that imply m = 3 are no longer present for inertial flows, see (Mei and Auriault 1991) article for more details.
In both the laminar and Stokes regions, the flow is steady. With increasing Re z the flow becomes unsteady defining the lower bound of the Transition region. Here, largescale oscillations are induced ) and ∫ | d dt |dV ≠ 0 where is any flow field variable and is the flow regime. The smallest spatial scales of the unsteady structures are, however, still comparable in size to the geometry preventing an isotropic turbulent energy cascade from taking place. In practice, this is visible as a turbulent energy spectrum decay curve with an exponent smaller than −5∕3 , as reported by Seguin et al. (1998).
In the Turbulent region, the scale of the flow instabilities is small enough to allow the existence of an isotropic turbulent energy cascade. Applying the Forchheimer equation, m goes up to 3 for a very short region and then returns to the second-order relation again; see (Lage et al. 1997) for details. In the current study, these regimes are considered with the main focus on the behaviour of the transition region.
The choice of Reynolds number affects the flow regime intervals; for Re z , the length-scale is the pore diameter, and the velocity is the mean streamwise velocity of the fluid, also called the interstitial velocity. A multitude of other examples can be found in Wood et al. extensive review (Wood et al. 2020). From now on we will refer to the particle Reynolds number Re p defined in Eq. (3) where U int is the streamwise interstitial velocity, D p is the particle diameter and is the kinematic viscosity.
For experimental work on macroscopic properties, the naturally occurring porous geometry is rather easy to represent. Usually, the porous bed can be chosen to be identical to the systems common in nature. Pore-scale flow visualization by PIV necessitates the usage of transparent bed particles and a spatial scale that the system is capable of resolving. Furthermore, a refractive index matching between the porous particles and the liquid is necessary for two reasons depending on what type of PIV-system is utilized. For 2D-PIV systems, the laser sheet is scattered out of the plane by refractive index variations. For tomographic PIV-systems, the backscattered light from the particles is refracted by the cylinders which cause inaccurate reconstruction of the velocity-field. This puts limitations on the type of porous materials and liquids that can be used. An extensive review of different possible refractively matched solid-fluid combinations are available in Wright et al. (2017).

Ordered thin porous media
Thin porous media as defined by Fabricius et al. (2016) can be classified into three types, mathematically expressed if h is the height of the cell and is the width then << h for HTPM (Homogenously Thin Porous Media). For PTPM (Proportionally Thin Porous Media) ≈ h and for VTPM (Very Thin Porous Media) >> h . In Fig. 1, HTPM, PTPM and VTPM are presented from left to right.
Ordered arrays of cylinders confined by two walls are of interest to study for at least two main reasons; the simplicity of the geometry makes it well-suited as a model problem, and 2D-versions of HTPM packings have been studied extensively by among others (Koch and Ladd 1997) and (Lundstrom et al. 2010). Secondly, flow between systems of tube bundles is common in several processes in nature and industry including heat exchangers (Khan et al. 2006), and flow through vegetation (Nepf 1999). Packed beds of spheres as studied in, for  can also be of interest for future studies of more complex geometries. It will, however, turn out that the flow fields in the geometries here studied become complex although the geometries are simple. The type of thin porous media that is investigated is PTPM since it has fully developed 3D-flow for all Re p . Two packings are considered, simple cubic and staggered arrays of mono-radii cylinders, see Fig. 2 for the elementary cells.
Conversion from the particle Reynolds number to the pore Reynolds number is possible by the relation from (Niven 2002) as which can be useful when comparing results in the literature. Here, is the void fraction of the bed, V p is the cylinder volume and A p is the cylinder area. There are multiple ways to choose the pore length scale d z , if the procedure above of Niven is utilized then Re p ≈ Re z for the cases in this study.

Earlier work
Examples of earlier experiments imaging flow fields in porous beds include planar PIV measurement through a randomly packed bed of spheres for laminar and turbulent flow, . Khayamyan et al. investigated the same packed sphere geometry for both laminar, and turbulent flow using both normal planar PIV  and stereoscopic PIV . For low Reynolds number flows, PIV has also been used to investigate flows through rectangular arrays of circular rods Agelinchaab et al. (2006). For flow on small scales, confocal microscopy has been applied to disclose 2D fluid velocities for flow through a porous medium of a sintered disordered packing of hydrophilic glass beads Datta et al. (2013). A similar porous medium was studied in Sen et al. (2012) with PIV yielding 2D low Re velocity fields. The flow characteristics of similar geometries have been studied numerically and mathematically by among others, Koch and Ladd (1997) The current study is an extensive continuation of the work presented by Larsson et al. (2018) where it was revealed that the time-averaged flow field changes considerably as Re increases. The differences compared to (Larsson et al. 2018) are the dimensions of the experimental cell, the bed porosity and the cylinder configuration. The data is also evaluated to a much larger extent revealing new insights.

Method
In this work, the experimental cell (made of PMMA plastic) has the internal dimensions 600 × 125 × 21 mm 3 , see Fig. 3. It contains flow guiding vanes, a mesh with 1.6mm square openings and a thread thickness of 0.5 mm, and a honeycomb with 2.5-mm circular holes to straighten the flow and promote uniform inlet flow. The uniform The quarter-cylinders at the corners of the staggered cell also have the radius r c inlet flow was confirmed using LDV measurements not described here. A valve, being located close to the outlet of the experimental cell and connected to a tube, is used to remove air bubbles entrained in the working fluid. This is necessary since bubbles in the flow influence the measurements and the results.
The porous medium is modelled with cylindrical rods made of quartz glass with a diameter of 15 mm and a length of 35 mm. The rods are placed on a plastic board with pre-drilled holes to keep them rigid and steady when fluid is flowing through the experimental cell. The plastic board is interchangeable meaning that different cylinder configurations can be used. The rods extend all the way to the upper wall of the experimental cell leading to a flow confined between two parallel plates. The bed porosity of the square arrangement is 0.6 implying that the side length of a unit cell is 21 mm, for the staggered cell (see right side of Fig. 3) the porosity is ≈ 0.63. A total of 60 (+1) rods (6 along the width and 10 along the length of the experimental cell) are used for the simple cubic arrangement, while 55 (+1) rods are used for the staggered arrangement. A notable difference from an ideal staggered packing is that the edge cylinders are missing, i.e. an optimal staggered packing would have half cylinders at the edges as well. This is known to affect the macroscopic porosity (i.e. the total void space) of the experimental cell. The rest of the description of the method follows (Larsson et al. 2018) and is to a large extent repeated here since the method is essential for the new results presented.

Refractive index matching
To be able to accurately measure through the quartz glass cylinders without optical distortion, it is necessary to match the refractive index of the working fluid to that of the quartz cylinders. The index matching was done following the procedure described in Larsson et al. (2018), with a mixture of mineral oil and heptane as the working fluid. The dynamic viscosity ( ) of the working fluids was measured with a rheometer (Bohlin CVO, Malvern Instruments & 2000 M/ ME rolling-ball viscometer, Lovis) and the density ( ) was measured through a Coriolis mass flow meter (CoriolisMaster FCB450, ABB), hence the kinematic viscosity could be determined. A pump was driving the flow and a cooling system in the tank kept the fluid temperature constant at In the measurements, fluorescent particles, PMMA spheres dyed with rhodamine B from MicroParticles GmbH (diameter between 20-50 μm , density of 1.19 g/cm3), were used as seeding. The particles emitted 560 nm light when excited with a 532 nm beam.
The refractive index matching procedure differs slightly since a cylinder is placed inside the experimental cell as opposed to outside in a separate cuvette as in Larsson et al. (2018), the refractive index matching can be monitored and fine-tuned at real experimental conditions. The fluid properties are summarized in Table 1. Two different mixtures with differing densities and viscosity but with the same refractive index at 25 • C are presented, mixture 1 is utilized for the tomographic-PIV experiments, while mixture 2 is used for the measurements of permeability. The density and viscosity are different between the mixtures since the mineral oil properties varied between the different batches used.

Data acquisition
The details of the tomographic PIV measurement setup is explained in detail in the work of Larsson et al. (2018). The system is commercially available from Lavision GmbH. It consists of a double-pulsed Nd-YAG laser from Litron (532 nm, 15 Hz, 200 mJ) with a laser guiding arm, four sCMOS cameras (5.5 MP, 16 bit, 6.5 x 6.5 μm pixel size), a programmable timing unit (PTU X) for triggering and synchronization, and a computer with the software Davis 8.4 to control the system and store the data. Volume optics generated the illumination. Figure 4 is an illustration of the experimental setup. The flow rates were varied to yield an Re p in the range of ≈ 10 to ≈ 800 for both arrangements. The measurements were carried out over a two-day period, first for the cubic cell and then for the staggered one. The flow-rates were successively increased from a value corresponding to Re p ≈ 10 up to Re p ≈ 800 . 250 image pairs were acquired for each Re p and the time between exposures was adjusted between 300 and 15,000 μs depending on flow rate. In addition to the PIV measurements, the pressure difference over the cylinder array was monitored with pressure transducers from General Electric (GE Druck Unik 5000). Pressure sweeps were performed by varying the frequency of the pump and continuously measuring the pressure drop across the cell.

Data processing
To ensure good reconstruction quality of the particles, the acquired images were preprocessed in several steps. A filter subtracting a sliding minimum (3x3 pixel window) was applied to remove background noise. The images were normalized with the local average to compensate for intensity differences of the laser illumination, smoothed with a 3 × 3 Gaussian filter and finally sharpened.
The calibration procedure described earlier generated a mapping function (based on a third-order polynomial function) between image and physical volume coordinates, and volume self-calibration, Wieneke (2008) was applied to correct and improve the mapping. Multiple iterations were performed to improve the fit of the mapping function and the final calibration error was below 0.08 pixels.
A fastMART (multiplicative algebraic reconstruction tomography) algorithm was used (6 iterations) in the reconstruction of the particle volume distribution, the illuminated volume was discretized with 2576× 2170×630 voxels. The velocity vectors were calculated by iterative, multi-pass FFT volume cross-correlation with a final interrogation volume size of 32 × 32 × 32 voxels with 75% overlap. Spurious vectors were removed and replaced, the vector field smoothed (Gaussian, 3 × 3 × 3 voxel) and finally the areas in the resulting images containing no information were masked out.

Mean convergence
To decide how many image pairs that are necessary to form a well converged mean value velocity-field, a convergence study has been carried out. There are two sources of variation between the snapshots, the first is the reconstruction error, which includes measurement errors originating from for example seeding density and camera noise (Sciacchitano 2019). The second source of variation is the true flow field changes due to a transient flow field. The reconstructions are averaged together forming a mean quantity and a standard deviation quantity for each velocity vector in the volume.
The residual of the convergence as a function of the amount of image pairs is presented in Fig. 5. The convergence residual Res is defined as where is the averaged quantity and k is the amount of image pairs. For most cases the residual is less than 10 −3 implying that the change in velocity between each residual is less than 0.1% of the mean. The value is closer to ≈ 1% for the staggered packing and Re p = 480 − 840 due to the central point being in the middle of a recirculation zone. From this, it is clear that 150 image pairs are enough to produce a wellconverged mean value. All results reported are therefore based on data evaluated from 150 consecutive image pairs.

Pressure measurements
Measurements of the pressure across the cell were carried out on the oil-heptane mixture 2, see Table 1. The points where the pressure is monitored is marked in Fig. 6. The sweep started at the maximum flow-rate and was decreased constantly to a minimum value while saving pressure data. Average values for each flow-rate were extracted and are presented in the following section.

Results and discussion
In this section, results from both the staggered and cubic arrangement are presented and a general overview of the flow features for different values of Re p is given by heat-maps and streamline visualizations. Also, key differences between the packings will be pointed out and discussed as the results are introduced.

Standard deviations of the flow field
The normalized standard deviation of the velocity is one measurement of how unsteady the flow is, see Table 2 where the average standard deviations of the flow field are normalized with respect to the average velocity. At an Re p of approximately 250 for the staggered and 380 for the cubic array the standard deviation values start to increase, indicating the beginning of the transition region where u ′ x ≠ 0. Here, std(x) is the standard deviation operator of x in the temporal dimension and < x > denotes a spatial mean of x taken across the whole measurement region. The term U int is the interstitial velocity as defined in the introduction.
As a result of the alignment of the laser light volume the region closest to the outlet/furthest downstream was not as illuminated as the centre and the region closest to the inlet/ furthest upstream. The data points were enough to form a well converged mean but the standard deviations are larger in this region compared to the surroundings, therefore this region is omitted from the standard deviation analysis. A convergence analysis is presented in Fig. 7, the residual is below 10 −3 for all cases.

Cubic array
At Re p ≈ 10 the flow approximates a Stokes-flow but recirculation zones start to appear behind the cylinders, see Fig. 8 where the central streamlines in the z-direction are coloured by the y-velocity and scaled by the interstitial Fig. 6 Points where the pressure was monitored, green is the higher pressure upstream of the cell and red is the lower pressure downstream Table 2 Average standard deviations of the flow field velocity scaled by interstitial velocity, indicates flow unsteadiness velocity, see Fig. 3 for coordinate system reference. The flow forms a central channel as Re p > 10 and shows no significant flow structure changes in the averaged field for the region 170 < Re p < 700 . However, at Re p ≈ 840 the lateral flow direction alternates. This is visible as yellow coloured streamlines at the first and third columns and purple streamlines at the second and fourth.
This alternating main channel flow/transversal motion was not observed in the less dense packing studied by Larsson et al. (2018). The impact of this effect on the macroscopic property of dispersion is expected to be significant. For numerical studies it has been standard practice to model single unit cells for the sake of computational efficiency, including Lundstrom et al. (2010), Koch and Ladd (1997) and many more. Since this effect alternates spatially in the   . 8 Visualization of the time-averaged central streamlines in the xy-plane, the streamlines are coloured by interstitial velocity scaled y-velocity. The cylinders are marked in black and the data location in the global cell is indicated by the green plane stream-wise direction at least two adjacent cells are necessary to represent it, therefore this result indicates that such a simplification may yield inaccurate results for some packings.
A partial explanation of the bidirectional flow may be provided by the phenomena of wall-jet attachment, this is called the Coanda effect and has been investigated by Vít and Maršík (2004) among others. In short, the wall-jet attachment increases with an increasing Reynolds number. Once the jet has partially attached the incompressibility of the fluid causes it to be pushed down at a velocity lateral to the main direction through the whole channel, this can be seen as the streamlines diverging as they hit the cylinders at Re p ≈ 840 in Fig. 8. The x-velocity distribution along the central-plane in the y-direction is presented in Fig. 9. As Re p increases the inertial core becomes increasingly flattened when the recirculation zones develop. At Re p ≈ 380 the re-circulation zones can be observed as symmetric dark lines in between the inertial core. At Re p ≈ 520 the re-circulation zones are visibly skewed, at Re p ≈ 840 the re-circulation zones are so skewed that the minima are no longer clearly observable. The dog-bone shape of the inertial core are in agreement with that presented in Larsson et al. (2018), although the packing density of the cylinders are higher in the current case.

Staggered array
At Re p ≈ 10 the flow in the staggered array approximates a stokes flow closely similar to the cubic configuration, see Fig. 10. As Re p increases to ≈ 360 the flow becomes chaotic and all recirculation zones are partially skewed or deformed. This phenomena was earlier observed by  at a similar Re p . This is also in agreement with the jump in the oscillations from a low value of approximately 10% equal to the background noise up to around 12.7%, see Table 2. It is worth noting that at these Re p a similar effect is not present for the cubic array. As Re p reaches the maximum value the average flow field once again becomes more periodic.
When scrutinizing the velocity distribution along the centre plane of the cell in the y-direction an interesting velocity field appears, see Fig. 11. From Re p ≈ 10 to Re p ≈ 80 the inertial effects are accentuated with the increase in Reynolds number and at Re p ≈ 170 the velocity field at the edges folds into a pair of epsilons with symmetry axes at the y-centre. The details of how this effect is generated are still unknown but the symmetry indicates that this effect originates at the walls and propagates towards the centre of the cell. This is confirmed by observing the previous and next centre-line in the stream-wise direction, see Fig. 12. The upper figure is one cell downstream and the lower is one cell upstream, in the downstream the inertial core has been destabilized/ skewed indicating that the effect is propagating in towards the centre of the cell. For Re p ≈ 170 , the average standard deviation of the x-velocity within this region is of the order 8.8% of the average x-direction stream velocity; this indicates that the effect is not the result of transient variation but rather a steady convective-laminar effect. As can be observed in Table 2 the Re p ≈ 170 flow has the same standard deviation fraction scales as the Re p ≈ 10 and Re p ≈ 80 cases, indicating a steady state. For physical beds imperfections are to be expected, this result indicates that such imperfections can cause flow field changes far downstream from the geometric disturbance. At Re p ≈ 250 the inertial core can be seen to start self-destabilizing, which is probably induced by instabilities generated in the lower velocity central region of the inertial core; at Re p ≈ 250 this region is visible as four blue dots in the two central inertial cores. For Re p > 250 the averaged structures in the flow remains prominent until Re p ≈ 790 where turbulence smooths out any laminar-steady/laminar-unsteady structures.
Another way to view the unsteady onset is by iso-surfaces of the velocity magnitude; the structure furthest to the right and left in the sub-figures of Fig. 13 marked by red rings corresponds to the velocity phenomena in Fig. 11. The variation of this structure with the Reynolds number can be observed in Fig. 13 where Re p from top to bottom is 170, 250, 360 and 790, the averaged structures in the channels are broken up by the unsteady fluctuations as the flow becomes turbulent. This is what causes the epsilon-effect to disappear with increasing Reynolds number in Fig. 11.

Pressure measurements
The Reynolds number Re p was varied between ∼ 80 to ∼ 1000 and 51 measurements of the pressure were carried out for the cubic array and 27 for the staggered array; see Fig. 14, the bars are scaled with the standard deviation. The pressure is scaled to be dimensionless by using eq (1); here k 0 = μL R 2 where is the dynamic viscosity of the liquid, R is the cylinder diameter, p is the pressure drop across the cell, see Fig. 6, and L is the experimental cell length. Variations of p * for this scaling of the pressure indicates deviation from Darcy's law.
The pressure curves approximately follow Darcy's law up until Re p ≈ 150 where the second order relation begins, again showing that from an engineering point of view Darcy's law can be valid to relatively high Re p , although the flow is not Stokesian, see Fig. 15.
For the cubic array the pressure is discontinuous at Re p ≈ 850 which indicates that some significant change in the flow field takes place, this is observed to coincide with the relative standard deviations of the u x velocity presented in Table 2 increasing sharply. It is reasonable to expect that the effects coincides with breakup of the inertial core by turbulent interaction. In order to scrutinize if this implies an alternating lateral flow velocity (as observed in Fig. 8), the velocity field at Re p ≈ 700 is visualized by streamlines in Fig. 16.
The recirculation zones are observed to be skewed and alternating for the first, rows but not as pronounced as for the  Re p ≈ 840 streamlines in Fig. 8. This indicates that the alternating effect may play a role in the pressure jump. It may also be noted that no jump in pressure can be observed for the staggered arrangement of cylinders in Fig. 14 in agreement with the more gradual transition into the turbulent region in Table 2. This further strengthens the hypothesis of inertial core skewing being the cause of this pressure jump.

Summary and conclusions
Tomographic-PIV measurements of the transition region were carried out on both a staggered and cubic arrangement of mono-radii cylinders. Pressure measurements of both cells were also carried out and the relation between pressure variations and flow features for the cubic arrangement of cylinders were discussed.
Details regarding how the pore scale effects observed such as the epsilon-effect in Fig. 11 or the alternating flow in Fig. 8 affects macroscopic properties such as the pressure drop or the dispersion are of interest for a continued study. Since the effect was observed to originate at the wall it also raises questions on how packing asymmetries may affect the overall flow features.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Fig. 15
Visualization of the time-averaged central streamlines in the xy-plane at Re p ≈ 170 with the cubic array to the left and staggered to the right, the streamlines are coloured by interstitial velocity scaled y-velocity. The cylinders are marked in black and the data location in the global cell is indicated by the green plane Fig. 16 Visualization of the central streamlines in the xy-plane for Re p ≈ 700 , the streamlines are coloured by interstitial velocity scaled y-velocity. The cylinders are marked in black and the data location in the global cell can be seen in Fig. 8